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ABSTRACT 


The role of forces produced by the musculotendon units in the stress develop- 
ment of the long bones during gait has not been fully analyzed. It is well known 
that the musculotendons act as actuators producing the joint torques which drive 
the body. Although the joint torques required to perform certain motor tasks can 
be recovered throu gh a kinematic analysis, it remains a difficult problem to deter- 
mine the actual forces produced by each muscle that resulted in these torques. As 
a consequence, few studies have focused on the role of individual muscles in the 
development of stress in the bone. 

This study takes a control theoretic approach to the problem. A seven-link, 
eight degrees of freedom model of the body is controlled by various muscle groups 
on each leg to simulate gait. The simulations incorporate Hill- type models of 
muscles with activation and contraction dynamics controlled through neural in- 
puts. This direct approach allows one to know the exact muscle forces exerted by 
each musculotendon throughout the gait cycle as well the joint torques and reac- 
tion forces at the ankle and knee. Stress and strain computed by finite element 
analysis on skeletal members will be related to these derived loading conditions. 
Thus the role of musculoskeletal dynamics and neuromuscular control in the stress 
development of the tibia during gait can be analyzed. 
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CHAPTER I 
INTRODUCTION 


The study of stress development in bone is an active area of research which 
is relevant to the prediction and prevention of injury. Stress fractures have been 
a recurring problem with military recruits during intense tr a i nin g and are now 
being seen more frequently in civilian athletes (Sharkey et al ., 1995). The need 
to better understand the development of these fractures is a driving force behind 
this research. While the problems being addressed by this research axe typically 
thought to reside in the fields of biomechanics or bioengineering, mathematics is 
central to the formulation and solution of these problems. In particular, control 
theory, fracture mechanics and numerical analysis are especially relevant in the 
approach developed in this investigation. 

Before issues such as failure of bone can be addressed, one needs to have a 
good understanding of the stress fields incurred in bone throughout tasks which 
are deemed normal. The integrity of the calculated stress fields during a task will 
ultimately depend not only on the accuracy of the model describing the material 
properties of bone but also on the proper interpretation and calculation of the 
loading conditions to which the bone is subjected. Defining the loading conditions 
which accurately reflect the state of a human bone during in vivo activity is a 
challenge due to the ethical limitations inherent in working with human subjects. 
The practicality of utilizing strain gauges bonded to the bone and other such 
monitoring devices is restricted due to the risk of infection and discomfort imposed 
on the subject. Thus non-evasive techniques are highly recommended. 

One of the goals of this research is to develop an analytical, predictive method 
for acquiring appropriate dynamic loading conditions for the long bones of humans 
throughout the gait cycle. Such gin analytical model would provide information 
that might otherwise be obtained only through evasive procedures. A promising 
approach to this problem lies in the forward analysis of a mathematical model 
capable of simulating normal gait which incorporates the interaction of the ac- 
tive and passive structures of the body. Muscles axe the active force producing 
components in the body while tendons, ligaments and bones are representative of 
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the passive structures. Although anatomic considerations imply the importance 
of considering the interaction between these structures, most studies in the past 
have concentrated on the behavior of these elements in isolation. Studies have 
shown that ligaments and bones behave differently in isolation by storing more 
energy, requiring more force to rupture, and sustaining increased elongation as the 
speed of loading is increased (Wainwright et al ., 1976). It is also hypothesized that 
stress fractures in bone represent either a failure of fatigued muscles to absorb im- 
pact (Markey, 1989) or are the result of uncoordinated muscular action (Lanyon, 
1984). Thus a reliable mathematical model must incorporate the combined ef- 
fects of the muscle- tendon-ligament-bone complex by considering both active and 
passive structures as biological control systems. 

Deriving loading conditions from a simula tion is a relatively new procedure 
and a chall enging modeling problem. Chapter II will emphasize advantages of this 
method. Chapters III and IV are devoted to the modeling of musculotendons. A 
great deal of importance is placed on the proper modeling of the musculotendons 
for two reasons. First, since the musculotendons are known to actuate the body, 
any model of gait which hopes to mimic nature must include a fairly accurate 
description of muscle. Secondly, it is a contention of this research that the muscu- 
lotendons play a major role in defining appropriate loading conditions of the bone. 
Since the validity of the loading conditions relies on the ability of the model to 
simulate normal gait as accurately as possible, Chapters V and VI will concentrate 
on the development of the musculoskeletal model, the incorporation of the mus- 
culotendon actuators, and the implementation of the simulations. Thus a major 
portion of this dissertation is devoted to defining and validating a model through 
which appropriate loading conditions on can be derived. 

The last portion of this dissertation will involve a finite element analysis of the 
tibia using the derived loading conditions. The effects of musculature in the load- 
ing conditions will be emphasized. These computations demonstrate how loading 
conditions from a forward analysis may be implemented in a stress analysis. More 
importantly, the results will indicate that the incorporation of both active and 
passive structures of the body are essential in accurately predicting the state of 
stress in the long bones of the lower extremities. 



CHAPTER II 

JUSTIFICATION OF THE APPROACH 


2.1 Introduction 

To study the stress and strain distribution on the long bones, one must know 
the loading conditions. The relevant static loading conditions used in the past 
include joint reaction forces, joint torques, bending moments, and in some cases 
estimations of individual muscle forces. The question to be addressed in this section 
is what type of analysis provides a good methodology for studying the effects of 
different loading conditions associated with various forms of gait. 

Typically the body is modeled as a multi-link object in this sort of analysis. If 
the model has n degrees of freedom (dof), then the equations of motion governing 
this object can be written in the following simple form 

MO = T + V + G + E (2.1) 

where 0 is a n x 1 vector cont ainin g segmental acceleration, M is the n x n 
“mass matrix,” and T,V y G y E are nxl vectors representing the internal segmental 
torques, inertial, gravitational, and external forces acting on the body. Once these 
equations are derived, there are essentially two approaches in finding the loading 
conditions derived from a certain motor task: an inverse-dynamic approach or a 
direct-dynamic approach. 


2.2 Inverse-Dynamics 

Inverse-dynamics is am approach that proceeds in a backwards manner, taking 
observed motions and then deriving from these motions the joint torques respon- 
sible. Thus motion acts as the input in this method and joint torques are the 
output. It is convenient to write equation (2.1) in the following form 

T = (MO - V - G - E) 


to emphasize the dependence of T on the body trajectories. This approach re- 
quires experimental, kinematic data, i.e., ground reaction forces, mass and inertial 
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characteristics of segments, and observed motions. If these kinematic variables are 
ass umed to be known functions of time, the set of differential equations of motion 
becomes algebraic and joint torques are easily computed (King, 1984). 

Although joint torques are easily derived, it is a nontrivial matter to discern the 
distribution of force among the muscles from this data. As Davy and Audu (1987) 
stated, “There are typically more unknown forces than can be determined in the 
equipollence relations between resultant and individual member force (Crownin- 
shield, 1978; Penrod et al., 1974), so muscle forces can not be determined directly 
from mechanical relations alone.” This is referred to as the “redundancy problem.” 
Typically the muscle set is reduced by grouping muscles of similar function (Pa- 
triarco, 1981) and by using EMG activity as a guide in determining which muscles 
were used during the motor task (King, 1984; Biewener, 1992). If the problem is 
still indeter min ate, then an optimization scheme is usually utilized; however the 
selection of appropriate optimization criterion or the construction of proper cost 
functions remains disputed (King, 1984). Some cost functions which have been 
used in the past include the sum of muscle forces (Seireg and Arvikar, 1973; Pen- 
rod et al, 1974), the sum of muscle stresses or a related quantity (Crowninshield, 
1978; Crowninshield and Brand, 1981), and energy expenditure rate (Hardt, 1978; 
Patriarco et al., 1978). All the works cited above included no excitation nor con- 
traction dynamics of the muscles, and thus are referred to as static optimization 
after Hardt (1978). Depending on the form of the cost function, the optimiza- 
tion problem is solved via linear programming, gradient-restoration algorithms, or 
another appropriate method (Davy and Audu, 1987). 

There are severed shortco mings associated with static optimization. These 
methods neglect the role of muscle dynamics and assume that muscle actions 
at any instant are independent of action at all other points (Davy and Audu, 
1987). Neglecting muscle dynamics often results in muscle force histories which 
are discontinuous in time (Yamaguchi, 1990). Furthermore, the results from static 
optimization are highly dependent on the cost function selected which sometimes 
predict physiologically unrealistic results (Hardt, 1978; Seireg and Arvikar, 1975). 
Adding physiological constraint equations can help prevent some of these unrealis- 
tic results (Pierrynowski and Morrison, 1985). Another drawback to this method is 
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that the predictions are hard to validate, since invasive techniques such as inserting 
force transducers are not applicable in a human analysis. Regardless of the short- 
comings, static optimization remains a frequently used tool in the determination 
of muscle forces. 

In s ummar y although inverse-dynamics is a valuable tool, it is limited in several 
respects. First, the determination of muscle force distribution from joint torques 
remains difficult, and the only validation of such a distribution is correlation with 
EMG recordings. Secondly the results from inverse-dynamics are highly depen- 
dent on accurate determination of joint angles and the calculation of joint torques 
(Patriarco et al . , 1981), but kinematic measurements are not always able to dis- 
cern subtle movements which might be critical (Pandy and Berme, 1987). Another 
limitation associated with this method is that it is not predictive in nature, that 
is, one is limi ted to studying motion which is actually produced by the subjects 
being monitored. Additional drawbacks associated with inverse-dynamics can be 
found in Hardt (1978). 


2.3 Direct-Dynamics 

In a forward or direct-dynamical analysis the joint torques are the inputs and 
the body motion is the output. Thus equation (2.1) is rewritten in the following 
form to emphasize this relationship 

$ = M~ l [T + V + G + E). 

It is important to realize what produces these joint torques. Joint torques are the 
accumulation of interned body forces such as ligaments, joint constraints, and of 
course, muscle forces. Muscles are the actuator in this system. Thus the true input 
into the system is indeed neural input. This neural input drives the muscles, and 
then the torques are an accumulation of both the passive and active structures 
which stabilize the body and propel it forward. Therefore controls for each muscle 
are needed to drive the system and the “redundancy problem” is revisited. Some 
type of optimization te chni que must again be utilized (Yamaguchi, 1992; Davy and 
Audu, 1987; Hatze, 1976; He, 1988; Levine et al., 1983). 
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Although optimization techniques must still be utilized in this forward analy- 
sis, the optimization techniques employed are typically dynamic, meaning muscle 
dynamics are incorporated into the model. Muscle dynamics and the equations of 
motion governing the body are necessarily coupled. If a certain motion, like gait, is 
to be simulated, then the cost function usually involves both a tracking error term 
and a term influencing the distribution of muscle force, like energy consumption 
(Davy and Audu, 1987 ). Once controls are achieved, then the system of differen- 
tial equations for the body and the system of differential equations governing the 
dynamics of each muscle group can be integrated forward in time to obtain the 
motion trajectories. Thus in a sense a direct-dynamic analysis is self-validating in 
that the controls specified do indeed result in the observed motion. 

In this research, we choose to use a direct-dynamic approach for several rea- 
sons. When using a forward analysis, the loading conditions needed for a stress 
analysis of bone, i.e., joint torques, joint bending moments, joint reaction forces 
and individual muscle forces can be instantaneously derived as the model pro- 
gresses forward in time. Since any motion is possible in a forward analysis, we 
are capable of generatin g an infinite amount of data. In addition, abnormalities 
in muscle function can be incorporated into the model. Thus fatigue or improper 
innervation in the musculotendon complex can be studied and the effects of such 
conditions on the integrity of the skeletal system can be quantified. In essence, 
we believe that a forward analysis provides the means to study a broader scope of 
loading conditions. 



CHAPTER III 

MUSCULOTENDON COMPLEX 


As stated before, muscle is the active force producing organ in the body. Any 
forward analysis of gait relying on neural input as the control must incorporate a 
fairly accurate model of muscle into the simulation. Thus a working understanding 
of the basic function of muscle is fundamental to this research. Before describing 
a muscle model appropriate for gait simulation, some of the basic properties of the 
musculotendon complex should be discussed. For a more detailed description of 
muscle properties the reader is referred to McMahon’s book, Muscles, Reflexes, 
and Locomotion," or Zajac’s Review, u Muscle and Tendon : Properties, Model, 
Scaling, and Application to Biomechanics and Motor Control." 

3.1 Architecture of Musculotendon 

Muscles are connected to the bone through connective tissues called tendons 
and can be oriented to the tendon in a parallel or oblique fashion. These muscle 
types are referred to as parallel or pennate muscles respectively. See Figure 3.1 
for the configuration of pennate muscles. All muscles in this study are treated as 
pennate muscles where a, the pennation angle, quantifies the pennation for each 
muscle group. Although pennation angle varies depending on the state of contrac- 
tion, the volume of muscle remains constant, a fact confirmed by Swammerdam in 
the 1660’s (McMahon, 1984). The proximal attachment of the musculotendon unit 
to the bone is called the origin while the distal attachment is called the insertion. 
The length of the musculotendon which is dependent upon body configuration is 
defined as the length along the bones from origin to insertion. 

The organization of muscle can be studied at various levels. Muscle generally 
refers to a bundle of muscle fibers which act in parallel. A single muscle fiber has 
a banded appearance. These bands divide the fiber into a series of sacromeres, 
the smallest functional unit of a muscle (McMahon, 1984). This striated appear- 
ance is the result of the different refractive properties of the A and I bands in 
the sacromeres (see Figure 3.2). At a microscopic level the sacromeres are seen 
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Figure 3.1: Geometry of Pennated Muscle. When muscles contract, they maintain 
a constant volume, and become more pennated. [Yamaguchi] 


to be a collection of Myosin and F-actin myofilaments, referred to as thick and 
thin filament respectively. This microscopic architecture will be relevant when a 
theoretical explanation of force generation is described. 

All levels of muscle structure down to the sacromeres behave in a homogeneous 
fashion. Because of this homogeneity, each level can be viewed functionally as a 
scaled version of another level. Each sacromere has the same relaxed length, and so 
the length of the muscle unit, L m , is a measure of how many sacromeres constitute 
one muscle fiber. Since each sacromere produces the same amount of force F a 
rm rW simil ar circumstances and muscle fibers act in parallel, total muscle force 
F" 1 should be a scaled version of F t . The scale factor between the sacromere and 
the mus cle uni t is generally proportional to the mean cross-sectional area of that 
muscle unit, defined as weight in grams divided by length in centimeters (Zajac, 
1989). 


3.2 Functional Properties of the Musculotendon 
3.2.1 Force-Length Properties of Muscle 

Muscle force is easily measured at various lengths under isometric conditions 
to produce force-length relationships. The curve produced when muscle is not 
stimulated is referred to as the passive force-length curve, f P (L m ). When muscle 
is fully activated, the curve that results is called the tetanized curve and must 
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Figure 3.3: Isometric Force-Length Relation for Muscle: (A) Full activation, (B) 
Active force scales with activation but passive is force is unaffected. [Zajac] 


be the result of both passive and active contributions. The difference in these 
two curves is referred to as tbe active force-length relation, The length 

at which the maximum active muscle force, F a , is developed is called the optimal 
muscle length, L a . See Figure 3.3 for the qualitative nature of these curves. At 
less than full activation, the force-length dependence is obtained by scaling the 
fully activated fi curve (Winter, 1987; Zajac, 1989). The reason for the passive 
behavior is easily seen as just a consequence of the elastic material properties of 
muscle and is independent of activation. 

A theoretical explanation for the active fi relation is called the Sliding Filament 
Theory. This model was developed simultaneously by H. Huxley and A. Huxley 
(no relation) in the 1950’s (McMahon, 1984), and is based on the microscopic 
nature of muscle (refer to Figure 3.2). It was proposed that during contraction the 
thick and thin filaments in sacromeres slide past each other. Force is generated 
in the crossbridges that occur between actin and myosin during the overlap of 



11 


6 5 4 3 2 1 I 

* tit* I j 



( M m)^Q99 »j« (.67 *+— 039— H 

t====~~===l ’ 


1 2 


| i n HiiiH lti l | 3 


ft ggjgjj g^ 5 



Figure 3.4: Active Tension-Length Curve for Muscles. [McMahon] 


these filaments. Crossbridges can only occur in the presence of Ca + . The calcium 
concentration is a byproduct of the activation level which explains why the fi 
relation is scaled by activation. These crossbridges occur, rotate to move the 
filaments closer and then detach (McMahon, 1984). A. Huxley later proposed a 
mathematical model of the Sliding Filament Theory based on the probability of a 
crossbridge occurring. This theory properly explained the // curve (see Figure 3.4) 
and why force is only generated in a nominal region, .5L 0 < L m < 1.5 L a (Gordon 
et al., 1966). This model has also been shown to predict many of the other features 
of muscle, for instance the force-velocity relation. 

3.2.2 Force- Velocity Properties of Muscle 

Active muscle force is also dependent on muscle velocity. When a muscle ac- 
tively shortens (concentric contraction) , it produces less force than it would under 
isometric conditions. A.V. Hill (1938) was the first to quantify this result with an 




12 


empirical hyperbolic relationship, 

(F" 1 + a){v m + b) = (F a + a)b. (3.1) 

So maximal active force is achieved when the velocity of the muscle, v m , is zero, 
and the velocity (when L m = L 0 ) at which active force is zero is called maximu m 
speed of shortening, v a . Maximum speed of shortening is the limit to how fast 
a single unloaded crossbridge in the sacromere of a specific muscles can shorten 
(McMahon, 1984). The quickness of a crossbridge shortening determines whether 
muscles are classified as fast or slow. 

In contrast to shortening, when a muscle is actively lengthening it is able to 
deliver forces above isometric forces. This relationship is not an extension of Hill’s 
equation to the negative region, as might be expected. In fact, Katz discovered 
in 1939 that the slope of the force-velocity relation when muscle lengthens is six 
times steeper for slow-lengthening than for slow-shortening. He also showed that 
there is a threshold which limits the amount of tension muscle can withstand. This 
threshold is approximately 1.8F 0 . At tensions beyond this, the muscle is known 
to “give,” a phenomena know as yielding (McMahon, 1984). A representation of 
the normalized force-velocity relation is depicted in Figure 3.5. This curve is also 
thought to scale with activation (Zajac, 1989). 

It is clear now that active force generation is dependent on three factors: length 
of the muscle, velocity of the muscle and the activation level. There is evidence 
to suggest that total active force generation is best described by a force-length- 
velocity relationship which is usually quantified as the product of the force-length 
and force-velocity curves (Zajac, 1989; Winter, 1987). Therefore it is convenient 
to vis ualiz e force generation as a surface in three dimension, where each activation 
level determines a new surface (see Figure 3.6). 

3.2.3 Force- Length Properties of Tendon 

The tendon is composed of material both internal and external to the muscle 
(see Figure 3.1). The most relevant property of tendon with respect to gait studies 
is described by the stress-strain relationship depicted in Figure 3.7. This curve 
will be used in a later section to derive a generic force-length curve for tendon. 
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Figure 3.5: Force- Velocity Relation for Muscles (C) Full activation when L m = L™ . 
(D) The force-velocity curve is also thought to scale with activation. [Zajac] 



Figure 3.6: Three-dimensional surface plot showing the dependence of active force 
generation on both length and velocity of muscle. [Winter] 
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Nominal Stress-Strain Curve for Tendons. [Zajac] 


It is common to assume that the stress- strain property of tendon is homogeneous 
throughout internal and external portions. The tendon tangent of modulus of 
elasticity increases with strain in the “toe” region (e* < .02) and then becomes 
linpar at higher strains until the tendon fails. Tendon failure occurs around 10% 
strain and 100 MPa (Zajac, 1989). The relevance of the interaction of tendon and 
muscle will be described in the next section. 



CHAPTER IV 

MUSCULOTENDON DYNAMIC MODEL 


4.1 Introduction 

There are a couple of reasons why the muscle and tendon should be grouped 
together and treated as the musculotendon unit in gait analysis. Tendon and mus- 
cle act in series and when the tendon stretches an amount approaching L m , the 
force-length characteristics of the grouped actuator differ significantly from the 
muscle alone. Thus it is important to treat the muscle and tendon as a muscu- 
lotendon complex especially when the ratio L\JL 0 is large, as is the case with 
muscles around the ankle and knee (Zajac, 1989; Hoy et a/., 1990). Another rea- 
son they should be regarded as one unit is because they function together with the 
dy nami cs of the body. In fact, the muscle, tendon and body segments constitute 
a coupled, multiple-input multiple-output feedback system (Zajac, 1989). Figure 

4.1 describes the interaction between these three units. 

As seen from the figure, dynamics of the musculotendon complex can be viewed 
as two processes acting in series. Neural input drives the activation dynamics and 
outputs the muscle activation which is an internal state quantifying the ability 
of the cross-bridge structure of muscle to generate active force. Then muscle 
activation as well the length and velocity of the musculotendon are inputs into 
the musculotendon contraction dynamics, the process which outputs muscle force. 
These muscle forces then drive the body dynamics which effects the length and 
velocity of the musculotendon complex. Thus the dynamics of the whole system 
are highly coupled. 


4.2 Musculotendon Contraction Dynamics 
4.2.1 Scaling and Curves 

When including several muscles in a model, it is advantageous to develop curves 
describing the attributes of a generic muscle. This model can then be scaled with 
appropriate parameters to reflect the dynamics of a particular muscle. In doing 
so, only four general curves need be specified: force-length function for passive 
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Figure 4.1: Block Diagram of the Coupled Dynmical Systems (A) Musculotendon 
actuators work with the body segments to produce the forces that propel the body. 
(B) D ynami cs of the actuator consists of two uncoupled parts: activation dynamics 
and contraction dynamics. (C) Musculotendon contraction dynamics represent the 
interaction between muscle contraction dynamics and tendon compliance. 
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muscles, force-length function of active muscle, force-velocity function for active 
muscle and the force-length relation for tendons. 

Scale parameters needed for each musculotendon group: 

1. Maximal isometric active muscle force: F a , 

2. Optimal muscle length: L a , 

3. Pennation angle when L m = L a : a 0 , 

4. Tendon slack length: L\. 

All forces will be scaled by F 0 and all lengths will be scaled by L a . Maximum speed 
of short ening will be defined as, v 0 = L 0 /t c . This quantity is used to render the 
muscle specific force-velocity relation. This quantity t c scales time and is known 
to vary for fast and slow muscles; however, r c = .Is is standard value used for all 
muscles types (Zajac, 1989). Note then that v a is not a new scale parameter in 
our model. 

Force-Length Relation for Passive Muscles (f p (i m )) 

This function represents the spring like behavior that muscles exhibit when 
activation is zero. A natural cubic spline which was fitted to data reported on 
Scott Delp’s wedsite (www.kin.ucalgary.ca/isb/data/delp/delp_mus) is utilized in 
this model. Figure 4.2 illustrates the spline as well as the reported data. The 
source of this data is well referenced (Delp and Loan, 1995; Delp et al., 1990, 
Delp, 1990; Zajac, 1989). 

Force-Length Relation for Active Muscles ( fi(L m )) 

This function represents the isometric force development in muscles when ac- 
tivation is 1, i.e., full activation. Note that force is only developed when L m is in 
the nominal range: .5 < l m < 1.5., as the “Sliding Filament Theory” suggests. 
Again a natural cubic spline was used which was fitted to data reported by Delp 
(Delp et al., 1990). Figure 4.3 details the correspondence between the data and 

the splines. 

Force- Velocity Relation for Muscle (/ v _1 (u m )) 

This curve accounts for the fact that muscle cannot change its length instanta- 
neously and that the contractile component is damped by a visous effect. Although 
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Figure 4.4: Inverse Force- Velocity Relation for Muscles: v m vs. F 


Hill’s equation (3.1) is believed to model force development when a muscle is short- 
ening, we choose instead to use a natural cubic spline which was fit to data collected 
while the muscle lengthened and shortened (Delp et al., 1990). Due to the way 
the curve is utilized in the formulation of the dynamics, it is convenient to express 
this relation in terms of the inverse /“ 1 (F). The curve is depicted along with the 
data in Figure 4.4. The output of the inverse will be normalized with respect to 
v c , the mA-rimum speed of shorte ning . In particular where F denotes a normalized 
active muscle force, 

v m = v m /v° = f; l (F). 

Force-Length Relation for Tendon 

A generic force-length relationship for tendon is derived by a method discussed 
by Zajac. Define tendon strain by e t = (L‘ - L\)jL\ and normalized tendon force 
by F l = F l j F 0 . Zajac assumes a generic force-strain curve (F l vs. e l ) based on the 
following two assumptions (see Figure 4.5): 

1. A no min al stress- strain curve can be formulated that represents all tendon 
(refer to Figure 3.7). 
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Figure 4.5: Generic Force-Strain Relation for Tendons: d 1 — F l vs. e l [Zajac] 


2. The strain in a tendon when force in the tendon equals the maximal isometric 
muscle force is independent of the musculotendon unit. We refer to this strain 
as e‘, and once this values is specified it determines a nominal value of stress, 
<j*. Although the value of each of these quantities is disputed, we will take 
e* = 0.033 and = 32 MPa (Zajac, 1989). 

Once these assumptions are agreed upon, the force-strain relation is achieved by 
scaling stress by <r‘ in the stress-strain curve and noting that 

** _ f*/a pt 

< Fo/A 

where A is the mean cross-sectional area of the tendon. Thus by scaling this 
generic force-strain relationship by F a and L‘ , a force-length function is found for 
the tendon. 

We use the following function to describe the normalized force-strain relation- 
ship 

, tN ( .10377 (e 9u ‘ - l) 0 < e l < .01516 

= 1 V / . (4.1) 

v ’ \ 37.526e t - .26029 .01516 < e* < .1 

Figure 4.6 show the correspondence between this function and data reported by 
Delp. If the strain in tendon reach values beyond .1, the tendon is known to 
rupture (Zajac, 1989). Since such an extreme value of strain should not occur 
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Figure 4.6: Force-Strain Relation for Tendons: F* vs. e l — 




during normal locomotion, this part of the curve need not be included in our 
analysis. 

Also derived from equation (4.1) is absolute tendon stiffness which defined by, 
K t ( F *) = dF t /dL t . Observe that 


K t {F t ) 


dF* 

dF l 

de t 

dF* 

' de l 

‘ dU 

F q 

dF* 


Li’ 

de* 



' (f?) 91( F* + .10377) 0 < F* < .3086 
(§) 37.526 F l > .3086 


(4.2) 


4.2.2 Derivation of Contraction Dynamics of Musculotendon Unit 

The dy nami cs for the model of the musculotendon unit depicted in Figure 4.7 
can now be derived by utilizing the curves defined in the previous section. The 
model used is referred to as a Hill-type model and represents the properties of 
the musculotendon as idealized mechanical objects. This model has been shown 
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Figure 4.7: Hill Type Model of the Musculotendon. 


to incorporate enough complexity while remaining computationally practical, thus 
qualifying itself as a good model for use in gait simulations (Audu, 1985). Muscle 
mass is assumed to be negligible in this formulation even though there are some 
stability issues which arise if muscle mass is absent in the dynamics (He, 1988). 

The length of the musculotendon L mt is defined to be the distance from origin 
to insertion, as stated previously. Note that this quantity is a function of the 
segmental orientations. From Figure 4.7, it is clear L™ 1 also satisfies the following 
equation: 


L mt = L t + L m cos a. (4.3) 

Muscle is known to maintain a constant volume. In two dimensions, this implies 
L w is constant, and so 


L w = L m sin a = L a sin a 0 . (4-4) 

Fqe is the force developed by the contractile element (CE), and Fpe is the force 
developed by the passive element (PE). Thus the force in the muscle is given by 
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F 171 = F C e + F pe = F 0 a(t)ML m )f v (v m ) + F 0 f p (L m ). 
Due to a force balance, 

F t = F m cos a. 

The dynamics for the musculotendon are expressed by 

"* dP dL‘ mK , (Ft)t , 


(4.5) 

(4.6) 

(4.7) 


dt dL‘ dt 

Now an expression for U can be formulated by differentiating equation (4.3), 

L* = L™* — (L m cosa — L m smaa). (4.8) 

Since L w is constant, differentiation of equation (4.4) yields an expression for d, 

L m 


0 = L m sin a + L m cos a a 
Substitution of a into equation (4.8) gives 


a = 


tana 


(4.9) 


L 1 - L mt — L m ( cos a + sin a tan a) 

fm 

= L"* - -=^. (4.10) 

cos a 

Since v m = L m , equations (4.5) and (4.6) can be used to solve for L m . In particular 
we have 


F l = F™ cos a = (F 0 a(i)/i(i m )/v(u m ) + F 0 f p {L m )) cosa 


and so 


v m = — = / 7 1 ( 

V 0 \ 


(FV cosa) - F 0 / p (Z m ) 
Foa(t)fi(L m ) 


)■ 


L m = = i 


, gp/cos a )-F 0 } p (L'") \ 

V F„a(()/,(L”) ) ' 


( 4 . 11 ) 

(4.12) 

(4.13) 
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Thus the differential equation describing the contraction dynamics of the muscu- 
lotendon is 

ir‘ - r* ] 

cos Of \ F 0 a(t)fi(L m ) )\ 

where 



£ = 


v 0 = 10L 0 , 


L m = 


- L { ) 2 + (L^y 


L\ (l + ^ 

L‘= V 


Lq Lq 

0 


91 


< F* < .3086 


. ^(l + ^raT) .3086 < F 1 
L mt - V 


cos Ot = 


F 1 

Ft = y 

•Co 


4.3 Activation Dynamics 

In this model of the musculotendon complex, neural excitement, u(t), is re- 
garded as a control variable which varies between 0 and 1. Neural excitement is 
related to contraction dynamics by activation, a(t), which scales the active force- 
length and force-velocity curves. Activation varies between 0 and 1, but due to 
the formulation of the contraction dynamics it is necessary to specify a minim um 
activation which is not zero, i.e., 0 < ^ < a(t) < 1. This process through 
which neural input is transformed into activation is called activation dynamics and 
is known to be mediated through a calcium diffusion process (McMahon, 1984). 
Although there is some evidence to suggest that activation dynamics are not inde- 
pendent of contraction dynamics (Zajac, 1989), we adopt the common assumption 
of independence. 
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Both quantities, u(t) and a(<), can be related to experimental data, i.e. EMG 
recordings. In this case, u(t ) is related to rectified EMG while a(t) is related to 
filtered, rectified EMG (Zajac, 1989). This relationship between u(t) and a(t) can 
be represented by the following bilinear form 


da(t) 

dt 


+ 


~ 1 
•Tact 


(j3 + { 1 - /?)«(*))] - a(t) = — • u(t) (4.15) 

J To ct 


where 


0 </?< 1 . 


Note that activation is fastest with time constant when u(t) = 1 and slowest 

with time constant To^//? when u(t) — 0. We took = .012s (Zajac, 1989), and 

f3 = .1 (Pandy et al ., 1990). 

As will be discussed, the controls, for the muscles incorporated in this model 
were adapted from Yamaguchi’s dissertation. To simplify the optimization tech- 
nique u tiliz ed in creating these controls, he assumed u(t ) to be piece-wise constant. 
With such an assumption, it is possible to find a closed form solution to the dif- 
ferential equation (4.15). In particular, assume 


Then, 


u(t) = u 0 U <t <tf 


a(t ) = u a + ( a(t T) - u 0 )e-^ (1+!>Uo)(t - to) U<t<t f (4.16) 


where 


a(t { ) = lim a(f). 

This simplification was useful in the simulation of our model since the number of 
differential equations to be solved is large and thus represents additional computer 
time. Figure 4.8 shows how the activation dynamics responds to a step input. Thus 
activation lags behind neural input with activation occurring at a more rapid rate 
them deactivation. 
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Figure 4.8: Activation Dynamics: a(t) vs. t. This graph represents the response 
of the activation dynamics to a step input of neural activation, i.e., u(t) = X[o,H - 



CHAPTER V 

THE MUSCULOSKELETAL MODEL 


5.1 Introduction 

In order to generate the loading conditions relevant during gait, it was neces- 
sary to develop a model driven by musculotendon actuators which could simulate 
normal and pathological gait. Although many gait models have been built, few 
included the complexity which is needed for realistic dynamic simulations. It was 
decided that a model developed by Gary Yamaguchi (1989) represented a good 
initial choice. With this model, it could be determined if the additional loading 
conditions provided by accurate descriptions of muscle forces would be beneficial 
in characterizing stress development in bones. Thus the following description of 
the musculoskeletal model is adapted from the dissertation work of Yamaguchi. 

Before a discussing the actual model, it is important to realize the original 
goal of the research pursued by Yamaguchi. In his study the feasibility of using 
functional neuromuscular stimulation ( FNS) to enable a paraplegic to walk with 
normal appearance and speed was questioned. FNS is a process where electrical 
currents are artificially applied to nerve and muscle tissue in order to stimulate 
muscles. In Yamaguchi own words, “The goal of the study was to examine whether 
minimal sets of muscles could be used in order to generate approximately normal 
gait trajectories without requiring either high levels of force or unduly precise 
control of muscle activation” (Yamaguchi, 1990, pg v). In order to study the 
feasibility of FNS, it was necessary for Yamaguchi to develop a model capable of 
reproducing most of the known “determinants” of gait as classified by Saunders et 
al. in 1953 while maintaining a practical limitation on the degrees of freedom as 
dictated by the technology available at that time. 

5.2 Skeletal Model 

The following is a brief synopsis of the model utilized in our research and 
the reader is referred to Yamaguchi’s dissertation for more details. The model 
constrains seven rigid-body segments which represent the feet, shanks, thighs and 


27 



28 



Figure 5.1: 3-D, 8 DOF Model Showing Segment Angle Definitions: (a) The stance 
hip has two DOF, while all other joints are revolute, (b) Front view showing pelvis 
list. Stance angles are specified with respect to the inertial jrame, n, while the 
swing angles are respect to the titled trunk reference frame, d. [Yamaguchi] 


trunk to 8 degrees of freedom. All joints are assumed to be revolute having one 
degree of freedom with the exception of the stance hip which has two degrees of 
freedom. This allows the hip to ab/adduct, a condition which reduces the degree 
of coupling between the trunk and the swing leg (Yamaguchi, 1989). Figure 5.1 
shows the generalized coordinates which were used to describe the configuration of 
the body. Joint angles and g 3 are measured with respect to the horizontal 

or transverse plane and the rotation of these joints occur about an axes parallel 
to n 2 . Movement of the stance leg is confined to the sagittal plane, but due to the 
extra degree of freedom granted to the stance hip, the swing leg and trunk can 
also move in the frontal or coronal plane through pelvic list. Joint angle <74 tilts 
the trunk as well as the swing leg about an axis parallel to fti = d\. Joint angles 
q 5 through q$ are measured in the titled plane and these rotations occur about an 
axis which is parallel to d 2 . 
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Figure 5.2: Definition for Inertial and Segmental Parameters: (a) Stance Leg 
(sagittal plane) (b) Swing Leg (titled plane). 


This musculoskeletal model represents a normal male with mass totaling 76 
kilograms. Figure 5.2 shows the definitions for the segmental dimensions and the 
inertial parameters used to specify this model. Segmental dimensions and inertial 
parameters are listed in Table 5.1 (Yamaguchi 1989). 

A key element of the model is that only one-half (14%, approximately left- 
toe-off, to 62%, approximately left-foot-flat) of the gait cycle is simulated in this 
analysis. The complete gait cycle (see Figure 5.3) can be reconstructed under 
the assumption of bilateral symmetry. Bilateral symmetry is not always valid 
since many asymmetries in gait have been reported (Winter, 1979); however, the 


I 


Table 5.1: Segment dimensions and inertial parameters used in this model (see 
Figure 5.2). All starred items are given with respect to the segments center o 
mass. (Yamaguchi 1989) 
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14 % ^ simulMlBd stmp — — ► 62 % 




Figure 5.3: Simulated Gait Cycle With Respect to Total Gait Cycle. This model 
assumes bilateral symmetry to reproduce 96% of the cycle. [Yamaguchi] 


advantages of this assumption warrant its inclusion. This assumption simplifies 
the modeling in that the stance leg is always the stance leg and the swing leg is 
always the swing leg. As a result, muscles in the swing and stance leg can vary 
according to the function of that leg in the simulation. With this assumption, it 
is also valid to have the stance toe constrained to the ground which eliminates 
one more degree of freedom. Thus this simplification requires less muscles to be 
modeled and less degrees of freedom. 
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Figure 5.4: Passive Joint Moment for Ankle Joint. Positive angles and moments 
correspond to plantarflexion of the ankle. 


5.3 Passive Constraints at the Joints 

The range of each joint angle was limited to normal ranges through the use 
of ligamentous constraints. If the joint angle stays within a nominal range, then 
the effects of the passive structure are min i m al, but when the nominal range is 
exceeded the passive torques grow exponentially. The general form of the passive 
moments is given by 


Mpas,, = h e~ k *<•“*> - A* -c9 (5.1) 

where 6 is the joint angle measured in radians , 6 is the joint velocity measured in 
radian per second and is measured in Newton-meters. Note that 6 2 <9 <6\ 

represents a nominal range for that joint. The parameters kj , 6j, and c were taken 
from Davy and Audu (1985, 1987) and modified by Yamaguchi (1989). Included 
in this passive moment is a damping term, —c9. This is vital to the model since 
a damping term was not included in the musculotendon model. See Table 5.3 for 
a list of the passive parameters and definitions of the joint angles in terms of the 
generalized coordinates. Figures 5.4 through 5.6 depict the passive joint torques 
under static conditions, i.e., neglecting the damping term. 
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moment (N m) 



knee angle (deg) 


Figure 5.5: Passive Joint Moment for Knee Joint. Positive angles and moments 
correspond to flexion of the knee. 


moment (N m) 



hip angle (deg) 


Figure 5.6: Passive Joint Moment for Hip Joint. Positive angles and moments 
correspond to extension of the hip. 
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Table 5.2: Coefficients for Passive Joint Moments 


ankle 

Q stance = Ql ~ + 5934 

0 swing = 2.1642 - ( <Z8 - <tl) 

fex = 2.0 
fc 2 = 5.0 
A: 3 = 9.0 
A: 4 = 5.0 

ci = 9.43 

0! = 1.92 (20° dorsiflexion) 

0 2 = 1.047 (30° plantarflexion) 

knee 

^stance *?2 ~~ 

fci = 3.1 
k 2 = 5.9 
k 3 = 10.5 
fc 4 = 11.8 

Cl = 3.17 

0! = 0.00 (full extension) 

0 2 = -1.92 (110° flexion) 

hip 

@ stance ^3 ^5 

0 swing = 3.14 - (?5 + ft) 

= 2.6 
fc 2 = 5.8 
k 3 = 8.7 
fc 4 = 1.3 

Cl = 1.09 

0 X = 1.92 (110° flexion) 

0 2 = 0.1744 (10° flexion) 
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Figure 5.7: “Soft” Constraint on Swing Leg During Double Support [Yamaguchi] 


5.4 Soft Floor Constraints 

Although the stance toe is constrained throughout the simulation, it was nec- 
essary to incorporate additional constraints in order to prevent the stance heel and 
swing foot from penetrating the “ground” and to eliminate excess sliding of the 
swing foot during double support. These additional constraints are considered to 
be “soft” (Hemami, 1975). The vertical ground reaction forces which acted on the 
heels of both the stance and swing leg were modeled as highly-damped, stiff linear 

springs 


Fnormal 


0 zheel > 0 

-(1.5 x 10 s ) zheel - (lx 10 z )zheel zheel < 0 


(5.2) 


where zheel is the height of the heel above the ground. On the swing heel an 
additional frictional force applied in a direction parallel to the ground in the sagittal 
plane is used to prevent excess sliding (see Figure 5.7). This force is proportional 
to the normal force applied to the swing heel, 


friction 


0 zheelnaing ^ 0 

— .^Fnormalizheeltunng)} zheel wing < 9 


(5.3) 


Once flat-foot is achieved by the swing leg a counterclockwise torque is applied 
to prevent the foot from penetrating the ground. The torque is model as a damped, 
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torsional spring which resists plantarflexion of the foot. This torque is given by 

Jo zheelmting > 0 or 5 > 0 

~ | -5696J - 38.19£ zheel^ing < 0 and S < 0 

where 8 is the angle the bottom of the foot makes with the ground. 

These “soft” constraints allow the same model to be used during the single 
a.nd double support phases of gait. If the ground had been modeled as “hard 
constraint, then one would lose a degree of freedom on the swing leg when the toe 
touches the ground. This would mean two models would be needed to simulate 
this phase of gait and a switch between the two systems would occur at heel strike. 

5.5 Incorporating the Musculotendon Actuators 

5.5.1 Minimal Set of Muscles 

As stated earlier, one of the goals of Yamaguchi’s dissertation was to find a 
minimal muscle set needed to simulate normal gait. This m i n imal set of muscles 
(see Figure 5.8) as well as the controls for each musculotendon actuator was deter- 
mined by a dy nami c optimization technique which will be discussed later. This set 
was found to be in good qualitative agreement with record EMG activity during 
gait (Yamaguchi, 1989). We will use the same muscle sets as his program derived. 
As a result, ten musculotendon units are incorporated into our model, five on 
each leg (refer to Figure 5.8). On the stance leg, the relevant muscle groups are 
the Soleus, Gastrocnemius, Vasti, Gluteus Medius & Minimus, and the Iliopsoas. 
While the swing leg utilizes the Dorsiflexors, Hamstring, Vasti, Gluteus Medius 
and Minimus and the Iliopsoas. Thus seven different musculotendon groups need 
to be specified in this model. Table 5.3 shows the constituent muscles composing 
each musculotendon group, and Table 5.4 gives a list of the parameters which are 
used to distinguish the dynamics of each particular musculotendon unit. 

5.5.2 Origins and Insertions 

In order to incorporate these musculotendon into the dynamics, we must place 
the muscles geometrically on the body segment so that the length of the muscu- 
lotendon, L mt and the velocity of the musculotendon, t» mt , can be derived as a 
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Table 5.3: Musculotendon Groups and Their Constituent Muscles (Hoy et al., 
1989; Yamaguchi, 1989) 


Musculotendon Groups 

Constituent Muscles 

Iliopsoas 

Iliacus 

Psoas 

Gluteus Medius 

Gluteus Medius 

& Minimus 

Gluteus Minimus 

Hamstring 

Semitendinosus 
Semimembranosus 
Biceps Femoris, long head 

Vasti 

Vastus Lateralis 
Vastus Medialis 
Vastus Intermedius 

Gastrocnemius 

Gastrocnemius Lateralis 
Gastrocnemius Medialis 

Soleus 

Soleus 

Dorsiflexors 

Tibialia Anterior 
Extensor Digitorum Longus 
Extensor Hallucia Longus 
Peroneus Tertius 



Table 5.4: Parameters De finin g the Musculotendon Actuators (Hoy et a/., 1989; 
Yamaguchi,1989) . 



Maximum 

Optimal 

Tendon 

Pinnation 


Isometric 

Muscle Fiber 

Slack 

Angle 

Musculotendon Actuator 

Strength 

Length 

Length 



pm 

o 

K 


a 


(N) 

(m) 

(m) 

(deg) 

Iliopsoas 

1474 

0.1269 

0.0850 

7 

Gluteus Medius &: Minimus 

2686 

0.0760 

0.0355 

10.4 

Hamstring 

2340 

0.4065 

0.3850 

8.7 

Vasti 

6482 

0.1096 

0.2250 

4.5 

Gastrocnemius 

1423 

0.0482 

0.4250 

14.8 

Soleus 

3599 

0.0243 

0.2700 

25 

Dorsiflexors 

1400 

0.1009 

0.2250 

6.9 
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Sfwrc# L*g 


Swing Lmg 



Figure 5.8: Minimal Set of Ten Musculotendon Actuators as Determined by Yam- 
aguchi 


functions of the state variables, i.e. the generalized coordinate, ft. As stated ear- 
lier, the attachment of the musculotendon to bone is specified by the defining an 
origin, the proximal attachment, and an insertion, the distal attachment. Effective 
origins and effective insertions are specified when the straight path from the actual 
origin to actual insertion passes through bone or out of anatomical range during 
certain body configurations. Origins and insertions are specified with respect to 
coordinate systems which are directed along the bones and are fixed with respect 
to the foot, shank, thigh and pelvis. The origin and insertion points for the 7 mus- 
cle groups used in this analysis are given in Table 5.5. We used the same origins, 
insertions and pathways as utilized by Yamaguchi in his model which were defined 
according Brand et al. (1982) and Hoy et al. (1990). 
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Table 5.5: Origins and Insertions of the Musculotendon Groups. X, Y and Z axes 
correspond to pelvic, femoral, tibial and foot coordinates system (see Figure 5.9). 
The foot reference frames is the same as the tibial reference frame at anatomical 
position. Subscripts a and e refer to actual versus effective origins/insertions. 


Musculotendon 

Point 

X 

(m) 

Y 

(m) 

Z 

(m) 

Reference 

Frame 

Iliopsoas 

O a 

0.0075 

0.1350 

-0.0400 

pelvic 


o e 

0.0260 

0.0293 

-0.0042 

pelvic 


U 

-0.0180 

0.3351 

0.0116 

femoral 


Gluteus Medius O -0.0155 0.0785 0.0076 pelvic 

& Minimus / -0.0159 0.3873 0.0589 femoral 


Hamstring O a -0.0409 -0.0455 -0.0140 pelvic 

I -0.0170 0.3800 0.0073 tibial 

Vasti a 0.0106 0.2026 0.0205 femoral 

I a 0.0170 0.3930 -0.0006 tibial 

Gastrocnemius O a -0.0203 0.0071 -0.0073 femoral 

I a -0.0368 -0.0429 0.0028 foot 

Soleus Oa -0.0268 0.2467 0.0006 tibial 

I a -0.0365 -0.0428 0.0056 foot 


O a 

-0.0155 

0.2175 

0.0134 

tibial 

o e 

0.0259 

0.0117 

-0.0093 

tibial 

h 

0.1035 

-0.0520 

0.000 

foot 




i 

3 










Figure 5.9: Definitions for the Foot, Tibial, Femoral and Pelvic Frames 


5.5.3 Calculation of Musculotendon Length, Velocity and Moment Arm 

The total length and velocity of the musculotendon complex, which is needed as 
input into the musculotendon dynamics, can be derived for most muscles through 
vector addition. In this case, the length of the musculotendon is given by 

L M = |Q^| + |OWt| + |£Jt| (5.5) 

where O and / refer to origin and insertion, and the subscripts e and a refer to 
actual and effective coordinates. The velocity of the musculotendon is then just 
the time derivative of the length. 
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Musculotendon forces were incorporated into the body dynamics as joints 
torques. When a musculotendon spans a joint J, a joint torque is realized across 
that joint . The joint torque due to musculotendons which do not span the knee are 
computed using standard vector cross product methods. In this case, the moment 
A | acting on the proximal segment is define as follows 


A ? = ±F 



(5.6) 


where p is a vector from the joint to a point on the line of action of the musculo- 
tendon, and the line of action is defined by the unit vector 


oX 

\XX\’ 

the tension in the musculotendon complex, F *, is produced according to equation 
(4.14), and the symbol x represents a vector cross product (see Figure 5.10). The 
cross product, 

_ o£ 

represents the effective moment arm associated with the musculotendon at a certain 
body configuration. The sign in equation (5.6) depends on whether the musculo- 
tendon acts to extend or flex the joint it spans. 

This general method works well for all muscle which do not span the knee. 
The complication which arises for those musculotendons which do span the knee 
(Vasti, Hamstring, and Gastrocnemius) occurs because of the complexity of the 
knee joint. In short, as the knee flexion angle varies this produces both a change 
in the location of the knee joint center and in the location of the patella. Since 
this complexity is not accounted for in the simple segmental model formulated 
here, it is best to use alternative methods to the vector methods discussed above 
when defining effective moment arms and the lengths of these musculotendons. As 
an alternative, we defined the effective moment arms of the Vasti, Hamstring and 
Gastrocnemius according to curves which were formulated by Yamaguchi from a 
planar model of the knee developed by himself and Zajac (1989). These curves 
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Figure 5.10: Muscle Pathway and Effective Moment Arm [Yamaguchi] 


define the moment arms as functions of the knee flexion angle (see Figures 5.11 
through 5.13). Thus the joint moments produced by these musculotendons which 
span the knee are given by 


M = ±F* • m e {9f) (5.7) 

where F l is produced according to equation (4.14) and m e is the effective moment 
arm as a function of the knee flexion angle, Of. 

We then calculated the length of these musculotendons by integration of the 
moment arm (Wendt and Johnson, 1985; Hoy et al. t 1990; Yamaguchi, 1989). In 
this method, the relationship between the length of the musculotendon (L mt ), the 
effective moment arm (m e ) and the joint angle is given by 


v 


mt 


dL mt 

dt 


m e (0f) 


d0 f 

dt 


Thus three more differential equations are added to the system. 


(5.8) 
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Figure 5.11: EflFective Moment Arm of Vasti as a Function of Knee Flexion 



Figure 5.12: Effective Moment Arm of Hamstring as a Function of Knee Flexion 



Figure 5.13: Effective Moment Arm of Gastrocnemius as a Function of Knee Flex- 
ion 
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5.6 Controls for the Musculotendon Actuators 

As stated earlier, when a direct-dynamic analysis of gait driven by muscu- 
lotendon actuator is to be simulated, musculotendon controls must be derived. 
Developing these controls constitutes one of the most difficult aspects in forward 
gait analyses. Thus one of the advantages of mimicking Yamaguchi’s model was 
that his controls, with min or alterations, could initially be used to produce loading 
conditions needed in our analysis of the long bones. 

Some form of dynamic optimization is usually utilized in developing controls for 
a forward analysis. Yamaguchi produced his controls through a two-phase process: 
(1) coarse optimization of controls and (2) fine-tuning via trial- and error. The 
first phase used the dynamic programming of Bellman (Kirk, 1970; Larson and 
Casti, 1978) to determine the mini mal set of muscles needed in the simulations 
and a crude estimation of muscle activation. During this process, Yamaguchi’s 
model was simplified significantly to ease the computational cost associated with 
dynamic programming; this included discretizing the state space. 

Dynamic programming offered Yamaguchi a couple of advantages over other 
optimization schemes. This approach allows the control to be dynamically opti- 
mized over the entire time interval as opposed to being optimized in a quasi-static 
fashion (Yamaguchi, 1989). Another advantage is that dynamic programming does 
not require the linearization of the dynamic equations of motion, nor do bounded 
controls present a problem. 

The cost function used in the dynamic programming consists of a error tracking 
term so that deviation from the nominal trajectory are penalized and a term which 
seems to be related to muscle fatigue (Crowninshield and Brand, 1981) 

Ji(k) = ]T w x> i{xi - x, idei ) 2 + £ w Ut i • ( 5 - 9 ) 

In equation (5.9), xj is one of the 2n elements comprising the state variable, X = 
(<7i, 92) <73> #4, <75t ?6) <77) Qs> Qu 92) 93, ?4i 95) 97» 9 s) (refer to Figure 5.1) associated 

with the body, F{ and PCS Ai are the force and physiological cross-sectional area 
of muscle l, and m is the number of muscles considered originally in 
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Figure 5.14: Control u{t) (dashed lines) and activation a(t) (shaded curves) for 
the ten musculotendon actuators used in Yamaguchi’s dissertation [Yamaguchi] 


Yamaguchi’s research. The nominal gait trajectories were contained in x^dea and 
were specified according to data recorded by Winters (1987) and Inman et al. 
(1981). 

Once a mini mal set of muscles (see Figure 5.8) and crude activation controls 
were formulated, Yamaguchi ran full model simulations and refined the neural 
controls. The controls Yamaguchi derived are pictured in Figure 5.14. When we 
ran our simulations, we started with these controls and proceeded to fine-tune 
them to meet the specific needs of our slightly modified model. As noted by other 
researchers (Yamaguchi, 1990; Pandy and Berme, 1987), it is surprising that such 
a complex, coupled and dynamically unstable system can simulate normal gait 
with such crude and simple controls. 


CHAPTER VI 
DYNAMIC SIMULATIONS 


6.1 Numerical Implementation 

Until recently a limiting factor on the complexity added to gait models was the 
difficulty associated with deriving the equations of motion. Kane’s vector based 
method for formulating equations of motion simplified this task, but still deriva- 
tions done by hand could still take months to perform (Kane and Levinson, 1985). 
Luckily Kane’s method lends itself to computer implementation. Programs now 
exist which are proficient in algebraic manipulations and Kane’s method (Symba, 
by Nielan, 1986; Autolev, by Schaechter and Levinson, 1987). These programs 
are able to calculate the entire set of dynamic equations in algebraic form for 
open-chain linked-segment models, making it possible to utilize models with more 
segments and more degrees of freedoms. 

Autolev was our choice of program because it provides a step by step approach 
to Kane’s method which allows some insight into the nature of the equations de- 
rived. As stated in the user’s manual (Kane and Levinson, 1996), “Autolev 
was created expressly to facilitate analyses based either on Kane’s method or on 
Newton-Euler equations.” In addition, Autolev produces efficient programs in 
Fortran or C for the numerical solution of ordinary, nonlinear differential and or 
non-differential equations. Thus Autolev was considered to be a valuable tool 
with regards to this research. 

A brief description of how Autolev formulates equations of motion may be of 
value. As stated earlier, Autolev allows the user to perform Kane’s step-by-step 
method online. Kane’s method, also known as Lagrange’s form of D’Alembert’s 
principle (Ju and Mansour, 1988), is based on dynamical equations of the form 

F, + F; = 0 (r = X, 2 p) (6.1) 

where F r are the constrained generalized active forces and F* are the constrained 
generalized inertial forces for a system S possessing p degrees of freedom in a 
Newtonian reference frame N. The constrained generalized active forces are defined 
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as 


F r = '£v?‘R i (6.2) 

1=1 

where v is the number of particles (or segments) that form S, Pi (i = 1, . . . , r) are 
the particles (or segments), vf* (r = 1 , . . . , p) are the constrained partial velocities 
of Pi in the reference frame N, and Ri is the resultant of all contact forces (such 
as ground reaction forces) and distal forces (such as gravitational forces). The 
constrained generalized inertial forces are defined by 

P; = XX' ■ R- (6.3) 

t=l 

where R* is the inertial force for Pi in reference frame N; that is, 

R* = —midi (6.4) 

where m* is the mass of Pi and a* is the acceleration of Pi in reference frame N. 

Equation (6.1) is equivalent to Newton’s second law of motion (Kane and Levin- 
son, 1996). Indeed if Ri is the resultant of all contact forces and distance forces 
acting on a typical particle Pi of the system S', and a* is the acceleration of Pi in a 
Newtonian reference frame N, then in accordance to Newton’s law the equations 
of motion are given by, 


Ri — midi = 0 (i = 1, . . . , v) (6-5) 

where m* is the mass of P{ and v is the number of particles of S. Dot-multiplication 
of equation (6.5) with the partial velocities v f* of Pj in N and subsequent summa- 
tion yields 


+ * (-™*«») = 0, (r = 1, . . . ,p), (6.6) 

*=1 i= 1 

which is indeed equivalent to equation (6.1) by the definition stated in equations 
(6.2), (6.3) and (6.4). 
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When using Autolev, one follows a simple recipe to obtain the equations of 
motion: 

• Declare a Newtonian Reference frame N. 

• Declare bodies, frames, points, and particles, specifying their mass and iner- 
tial characteristics. 

• Declare generalized speeds and their time-derivatives. 

• Declare generalized coordinates and their time-derivatives. 

• Create the kinematic differential equations by relating the time-derivatives 
of the generalized coordinate to the generalized speeds. 

• Form position vectors and directional cosine matrices. 

• Form angular and linear velocities , which AUTOLEV can do for you. 

• Impose motion constraints, if necessary. 

• Form angular and linear accelerations , which Autolev can also do for you. 

• Enter expression for forces and torques. 

Once these components are entered, Autolev can generate the equations of mo- 
tion and Fortran or C code which will integrate these equation forward in time 
using a fourth-order Runge Kutta integration scheme. In our model, the influences 
of the muscles where realized as torques entered about the joints, and ground reac- 
tion forces were entered as forces applied to the heels. Most of the components of 
the musculotendon dynamics were entered into Autolev directly; however, some 
issues, such as the controls and the spline approximations to various curves, were 
entered directly into the C code in the form of subroutines. 

One of the additional benefits in using AUTOLEV was its ability to find instan- 
taneous loading conditions needed in our finite element code to analyze the stresses 
in long bone. Once the equation of motion were generated, one could specify ad- 
ditional generalized speed, properly constrain them and then solve for the contact 


Figure 6.1: 2-Dimensional Representation of Our Gait Simulation 


or joint reaction forces found at the distal and proximal ends of the segments. 
Autolev was also able to resolve the tangential components of the muscle forces 
along the bones, which enabled us to specify surface tractions along the bone. 

6.2 Dynamic Simulation 

Most simulations were performed on the Digital Alpha workstation. The av- 
erage time needed to perform a full simulation was approximately 3 hours. Many 
simulations were ran before the model performed appropriately. Figure 6.1 shows 
the graphics corresponding to one of the better simulations. Notice that the fig- 
ures do appear to mimic normal gait, with a few exceptions which will be discussed 
later. 

For obvious reasons, the simulation and many of the results reported in this 
section me very similar to the results reported by Yamaguchi. However, since our 
models do differ slightly, e.g., we used different muscle models, slight variations do 
occur and shall be pointed out to the reader. 

6.2.1 Controls and Initial Conditions 

As stated earlier, one of the reasons for utilizing Yamaguchi’s model was that 
the basic controls laws needed for a gait simulation were already formulated. We 
thought that this would accelerate the rate at which we could arrive at the loading 
conditions needed for our continuum analysis of the bones. In retrospect, it is 
not clear that this did indeed speed up the process. This is because it is hard 
to mimic such a complex process when all details are not clearly defined. The 
degree to which small variations in the model affect the overall performance of the 
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simulation is hard to define, and one is left wondering if variations are due to the 
slight differences or if a fundamental element has been neglected or misinterpreted. 

With that said, we found that our final model did indeed mimic the overall 
performance of Yamaguchi’s model quite adequately, and our control laws varied 
only slightly from the original set as published by Yamaguchi (1989). The final 
version of our controls for each musculotendon unit are illustrated in Figure 6.2. 
The reader is referred back to Figure 5.14 to evaluate the differences between these 
two control sets. 

One possible explanation for some of the differences could be the variation of 
the muscle model used. In Yamaguchi’s dissertation, there is a schematic diagram 
depicting the Hill type musculotendon model he utilized. In this diagram, a se- 
ries elastic element within the muscle is depicted; however, no further information 
on the approximate form of this elements was given. As a result, we choose to 
ignore this feature. Thus our musculotendon actuator undoubtedly reacts some- 
what differently from his. It should be noted that Zajac gives some fundamental 
reasons which validate the exclusion of this element, in that it’s inclusion raises 
issues about the constituency of the basic notion that sacromeres and fibers act 
in concert. Zajac also states that in most instances tendon compliance dominates 
and that the elastic element in muscle can be neglected (Zajac, 1989). Thus we 
felt justified in the exclusion of this element from our musculotendon model. 

Another possible explanation for the differences needed in our simulation could 
be the result of the variance in the initial conditions used to start the simulation. 
In order to start the simulations, the program must be provided with the initial 
segmental orientations (<ft, the initial angular speeds (<ji, the initial 

force in each musculotendon, and the initial lengths of the musculotendons which 
spanned the knee joints whose lengths were found using the moment arm integra- 
tion method. Since initial conditions were not specified in Yamaguchi’s disserta- 
tion, we took educated guesses at what the model’s initial segmental orientation 
should be by studying Yamaguchi’s figures and graphs and then by referring to the 
gait data reported by Winter (1987). The initial force in each musculotendon was 
found by running mock simulations to find the steady state muscle force achieved 
by the initial control when the motion of the model was constrained. Finally the 
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initial lengths for the musculotendons which crossed the knee were determined by 
trial and error. We sought an initial length of the musculotendon such that the 
muscle maintained an appropriate length throughout the gait cycle. The initial 
conditions defined for the angular velocity of the segments was crucial for a proper 
simulation of gait. Of particular importance was the initial sagittal plane velocity 
of the HAT. Due to its enormous mass, errors here were hard to overcome. Table 
6.1 defines our initial conditions. 

6.2.2 Joint Angles and Joint Torques 

A standard means of comparing and validating gait simulations is achieved by 
displaying the joint torques which drive the system and the resulting joint trajec- 
tories. Figure 6.3 depicts the standard definitions of the joint angles. Notice that 
these joint angles are distinct from the generalized coordinates used in the analy- 
sis. Figure 6.4 displays the resulting joint trajectories in our simulation. Although 
these trajectories appear to be in good qualitative agreement with the trajectories 
reported by other researchers (Yamaguchi, 1989; Winter, 1987; Inman et al., 1981), 
there are some critical areas of concern. One such area is the slight hyperextension 
of the stance leg knee joint. Another area of concern is the exaggerated flexion of 
the hip. 

Joint torques represent the combined effects of the active and passive structures 
which drive the body. In an effort to describe all the components which accumulate 
in the total joint torque, two graphs are displayed for each joint angle (see Figures 
6.5 through 6.9). The first graph displays the total joint torque and distributions 
of the active and passive components. The second graph breaks down the active 
component, i.e., the muscular component, into the individual contributions from 
each musculotendon unit sp ann ing that joint. The second graph thus gives one 
an idea of the role played by each musculotendon throughout the gait cycle. It is 
important to understand that these pictures do not give the whole story in that a 
torque produced by a muscle spanning one joint can and usually does accelerate 
all body segments to a lesser or greater extent. This idea of describing muscular 
function by determining the segmental accelerations caused by each musculotendon 
was suggested by Zajac and Gordon (1989). We find that our total torques are in 



Table 6.1: Initial Conditions for the Gait Simulation 
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Figure 6.3: Joint Angle Definitions 


good agreement with Yamaguchi’s reported joint torques, and these torques are in 
general accordance with the data reported by other researchers (Davy and Audy, 
1987; Hardt, 1978). 

6.2.3 Musculotendon Analysis 

In this section, the dynamics of each musculotendon unit will be reviewed in 
terms of the normalized force and power produced by each musculotendon. Figure 
6.12 reports the normalized force realized in the tendons of all ten musculotendon 
actuators. These normalized forces, F t /F 0 , were obtained through the dynamics 
described in Chapter V. Although Yamguchi did not report such results, similar 
forces are predicted by other authors (Pierrynowski and Morrison, 1985; Brand et 
a/., 1986). 

Another means of describing the functions of the muscle, tendon and the com- 
plete musculotendon unit is through the analysis of the power trajectories. These 
trajectories are depicted in Figure 6.13. This power reflects the rate at which the 
muscles and tendons are expending and absorbing energy. When a muscle or ten- 
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Figure 6.4: Joint Trajectories 
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don produces force while shortening, a concentric contraction, energy is released 
to the system. In contrast, if force is produced while the unit is lengthening, an 
eccentric contraction, then energy is absorbed or stored. This gives insight into 
the function of the musculotendons. 

To exemplify these ideas, we will examine the power histories of the plantarflex- 
ors, the Soleus and Gastrocnemius, of the stance leg. After heel strike of the stance 
leg and during single leg support the total power, the summation of the powers in 
the tendon and muscle, of the plantarflexors are predominately negative as they 
store energy or absorb the weight and maintain the upright position of the leg. 
However, during the push off phase this power becomes positive as they release 
this energy to facilitate lift off of the stance leg and the shift of weight to the swing 
leg during double support. In addition, the hamstring is seen to absorb energy as 
it brakes the extension of the swing leg. 


{Normalized Force} 


{ Soleus } { Dors i flexors} 




{Gastrocnemius } 



{Hamstring} 



{ Vasti } 



{ Iliospoas } {Iliopsoas } 




Figure 6.12: Normalized Force Histories for each Musculotendon. Left column 
refers to stance muscle while the right column refers to the swing side muscles. 
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Figure 6.13: Normalized Power Curves for each Musculotendon. Left colu mn refers 
to stance muscles while the right column refers to the swing side muscles. Positive 
power corresponds to the rate of energy the muscle expend while performing a 
concentric contraction. Negative power is the rate of energy the muscle absorbs 
during an eccentric contraction. All curve are normalized by the weight of the 
model, 76 kg. 


CHAPTER VII 

A CONTINUUM ANALYSIS OF SKELETAL ELEMENTS 
7.1 Introduction to Bone Mechanics 

Before giving the results of our continuum analysis of the long bones, it might be 
beneficial to give a brief summary of some of the issues involved in bone mechanics. 
First some of the difficulties associated with modeling the skeletal structures will 
be defined. A brief synopsis of the work which has been done in this field will then 
be presented. This summary will include the common simplifying assumptions 
under which bone is modeled. In conclusion, special attention will be given to the 
research which has attempted to include some of the effects of the musculature in 
the analysis. 

Bone is known to have a complex geometry and nonhomogeneous constitution 
and these factor are responsible for many of the difficulties associated with defining 
the mechanical properties of the skeletal members. One only needs to analyse 
the structure of the femoral head or pelvic girdle to appreciate the nonstandard 
geometry inherent in bone. With regards to the nonhomogeneity of bones, it 
is known that bone is the composition of three different materials. The hard 
outer shell is referred to as cortical or compact bone while the interior of the 
bone is composed of cancellous bone, also referred to as trabecular or spongeous 
bone, and bone marrow. Trabecular bone is a very complex material which is 
noncontinuous even on a macroscopic level and is best described as a structure not 
a material (Huskies and Chao, 1983). Thus bone is most accurately described as 
an anisotropic, nonhomogeneous material. 

Before giving a brief review of some of the current work, it should be noted 
that bone mechanics is not a new area of research. Some of the classical studies 
which are cited often are the products of Wolff (1870, 1892) and Koch (1917). 
Wolff postulated that the patterns seen in the trabecular of bone are aligned with 
the principal stress fields developed in the bone. This hypothesis became known as 
“Wolff’s Law” and research continues to be conducted in an effort to substantiate 
this idea (Hayes et al . , 1982; Koch, 1917; Rybicki et ai, 1972). Later, Koch did 
some ground breaking work as he gave the first detailed geometric description of 
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the femur and went on to calculate the stresses induced by loadings which were 
assumed to occur during normal activities. This work was a landmark in that 
it was the first attempt to model the bone as a beam and the first to attempt to 
include some effects of the musculature on the stresses and strains. Although Koch 
underestimated the joint and muscle forces (Duda et a/., 1997), he still concluded 
that that muscular action could reduce the bending stresses (Koch, 1917; Rybicki 
et at ., 1972). 

Since then many studies have followed Koch’s lead by modeling the bone as 
a beam. Viano et at. (1976) assumed that the femur was isotropic and modeled 
it as a hollow cylinder. While Piziali (1976) modeled the tibia as a cantilever 
beam fixed at the knee and loaded at the ankle. He suggested that bone be 
viewed as transversely isotropic, that is bone displays one set of elastic properties 
in one direction and a second set in the perpendicular plane. Thomsen (1990) 
modeled the tibia as straight, twisted non-uniform Timshenko beam. He viewed the 
compact bone and cancellous bone as homogeneous, linearly elastic, transversely 
isotropic material, while bone marrow was consider completely flexible. Thomsen 
concluded that the stiffness of the trabecular can be ignored, the twist of the 
tibia is insignificant and thus that a simplified approach modeling the tibia as a 
unif orm beam using mean cross-sectional properties can be justified as long as shear 
deformations are considered. There are many more authors who have modeled the 
long bones as beams (Cowin et al ., 1985; Huskies, 1983; Salathe et at ., 1989). 
This should give one a feel for the general assumptions made while studying the 
mechanical characteristics of bone. 

Finite element techniques were introduced into the field of bone mechanics in 
1972 about fifteen years after their conception into engineering mechanics (Huskies 
and Chao, 1983). The usefulness of this technique resides in its ability to han- 
dle complex geometries and nonhomogenous materials. For a nice review on the 
progress of FE techniques in orthopedic biomechanics throughout the first decade 
of its use, the reader is referred to the paper written by Huskies and Chao (1983). 
Rybicki et al (1972) did a comparison of beam theory and a continuum theory 
in the form of a finite element program when he studied the femur. This analysis 
is of particular interest to this dissertation because Rybicki did study the effects 
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of certain muscle groups on the distribution of stress and internal moments. He 
concluded that while classical beam theory is appropriate in the shaft of the femur 
it is inaccurate for regions of the greater trochanter, femoral head and areas of 
muscular attachment. He also stated that muscles have a pronounce effect on the 
maximum stress and total strain energy of the femur. Thus we are encourage to 
pursue a continuum analysis. 

As stated above, there is particular interest in the research that has attempted 
to included the effects of the muscular system into the loading criteria. This 
dissertation is interested in the response of bone to natural settings where muscular 
activity plays a role. Unfortunately, it is often hard to define these natural loading 
conditions. The determination of muscular forces is nontrivial, as was discussed in 
Chapter II. As a result, most studies have focused on the in vitro response of bone 
to artificially applied loads (Collier et al. } 1982; Curry, 1959; Viano et al., 1976; 
Reilly and Burstein, 1975); however, there are some studies which have attempted 
the inclusion of normal loading conditions (Duda, 1997; Koch, 1917; Lanyon and 
Smith, 1970; Orne and Manke, 1975; Toridis, 1969; Sharkey, 1995, Rybicki, 1972). 
The studies of Koch and Rybicki have already been mentioned, but it should be 
noted that the muscular forces they utilized were not derived from a particular 
task. Toridis (1969) developed a theory in which the bone is modeled as a three 
dimensional elastic beam where muscular force can be included along the bone as 
point attachments. Although the theory of how to compute these stress field is 
presented, a detailed example of its application is absence and the reader is not 
informed on the effects of the muscular forces into the analysis. Sharkey (1995) 
conducted an in situ experiment studying the metatarsal strain in cadaveric feet 
in an attempt to characterize the development of metatarsal stress fractures which 
occur frequently in the army. By simulating physiological loading and contraction 
of the plantar flexors, he was able to conclude that the fatigue of the flexors have 
a profound effect of the development of stress fractures. Lanyon and Smith (1970) 
conducted in vivo experiments by placing strain gauges on the tibial shafts of sheep. 
Among their conclusions is the statement that local deformation that culminate in 
fracture are in many case the result of uncoordinated muscular action. Fortunately 
or unfortunately, these types of studies are unethical in human analysis, and one 
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must resort to other means to derive the loading conditions for normal activity. 
Duda (1997) in his paper entitled, “Internal Forces and Moments in the Femur 
During Walkin g,” attempted to do just this. In this paper, he included the effects 
of “all thigh muscles, body weight and contact forces.” The muscle forces were 
taken from the work of Brand (1982, 1986). Thus were derived using a “static” 
optimization technique as described in Chapter II. Although Duda states that this 
was not a continuum analysis, it is unclear from the paper how he modeled the 
femur. Regardless, he does conclude that the muscles play a substantial role in 
balancing the loads within the femur. 

7.2 Finite Element Analysis of Tibia 

It is a contention of this research that the interrelationship between neural 
control, applied muscle forces and the resulting motion of body segments influences 
the development of stress in skeletal members. To explore this notion, a continuum 
analysis of the tibia will be carried out in which the boundary inputs are derived 
from the results of the dynamic model of gait. The main point of the approach 
that is advanced here is that direct dynamics provides instantaneous information 
regarding the direction and magnitude of muscular forces, net moments at joints 
and joint reaction forces. This information will be utilized in a stress analysis of 
the segments that represent the various bones. 

To illustrate the mathematical issues involved, a general formulation of the 
relevant boundary value problem is first presented. In particular, let Cl denote a 
body in R s with boundary T = T t U Tj, T t fl = 0. For u, x 6 R s , let u(x , t ) 
denote the displacement of points x 6 Cl at time t. The internal forces, more 
precisely the stress in Cl, may be described in terms of a tensor a that satisfies the 
equation of motion 


div(cr) + f(x, t) = pu u (x, t) x E Cl. (7.1) 

Here f(x, t) denotes externally applied forces. In the context of the biomechanical 
situation consider here, Cl represents bone and / corresponds to weight or exter- 
nal surface tractions which are derived from muscular forces. The displacement 
field must satisfy the equation of motion subject to boundary conditions that are 
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typically of the form, 


an — t, ier ( 
u — /i 0 , 16^ 

where n denotes an exterior unit vector to T t . Additionally, one must specify 
a constitutive relation that relates the states of stress to the deformation in fi. 
Though bone is most accurately modeled as anisotropic material, preliminary work 
has been carried out under the assumption that bone is isotropic. In this case the 
components of stress, < 7 y, i,j — 1, 2, 3 may be expressed in terms of displacements 
by ‘Hooke’s Law’ which states, 

&ij = 2 fieij 4 - Y~^ Sij ( €n 622 633 (^- 2 ) 

In the above equation, /x and 1 / are elastic parameters and the ‘strain’ is defined 
by 



Specific problems that correspond to simplified modes of deformation will now be 
presented. 

In the event that a bone experiences a pure tensile or compressive deformation 
and inertial effects are negligible, it suffices to solve the equation of equilibrium, 

^(E(x)A(x)^\+f(x)=0, 0<x<l (7.4) 

where E(x ) denotes Young’s modulus, A(x) the cross sectional area of the bone 
which has length /, /(x) is an applied surface traction, e(x) = ^ is the strain and 
a(x) = E(x)e(x) is the stress. This equation must be solved subject to boundary 
conditions that reflect the joint reaction forces and specified displacements at the 
joints. 

In this study, the swing leg tibia is analysis at heel strike and flat foot. At these 
two instances of gait, it is believed that the tibia can be described as undergoing 
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with surface tractions 



Figure 7.1: Axial Stress Induced During 
compression when muscular attachments 
muscular attachments are neglected, the 
in scales between (A) and (B). 


x 10° without surface tractions 



Flat Foot (A) Note that the tibia is in 
are included in the analysis. (B) When 
bone is in tension. Note that difference 


a pure tensile or compressive deformation. The results in Figures 7.1 and 7.2 
are finite element approximations to a boundary value problem that correspond 
to solving equation (7.4) in which f(x) represents the surface tractions due to 
muscular loading, in particular the tangential component of the vasti, hamstring 
and dorsiflexor muscle forces. It was assumed that the distal end of the tibia was 
fixed at heel strike and flat foot and the boundary condition at the knee joint was 
given by the normal component of the joint reaction force, f Q . In particular, the 
boundary conditions are 

du 

n(l) = 0 £(0).4(0)— (0) = U (7-5) 

A point to be emphasized is that the inclusion of muscular effects in boundary 
inputs significantly alter the distribution of stress throughout the bone. This con- 
clusion is supported by other researchers (Duda, 1997; Lanyon and Smith, 1975; 
Rybicki, 1972; Sharkley, 1995). With the inclusion of muscular attachments por- 
tions of the bone are in compression. This is in good qualitative and quantitative 
agreement with the results of Lanyon (1974) and Duda (1997). This should be 
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Figure 7.2: Axial Stress Induced During Heel Strike (A) Note that the tibia is in 
compression when muscular attachments are included in the analysis. (B) When 
muscular attachments are neglected, the bone is in tension. Note that difference 
in scales between (A) and (B). 


contrasted with the tensile stresses that are predicted when muscular tractions are 
neglected (see Figure 7.2). The results are not surprising in view of the fact that 
muscular forces often exceed body weight and thus dominate the loading on bones 
both at the joints and at the muscle attachment area (Hardt, 1978). 

Undoubtedly, more complex modes of deformation must be addressed in order 
to accurately represent the stress in the bone. Ultimately we intend to utilize 3-D 
finite element solvers to conduct the stress analysis of the skeletal components. 
However, a more immediate approach that incorporates the effects of axial, bend- 
ing and torsional loading will be based on a model of the skeletal system as a 
configuration of beams in a 3-D structure referred to as a frame. The work done 
by Toridis (1969) and Rybicki (1972) could serve as templates for this analysis. 
In this approach, the require boundary inputs would include the distribution of 
muscular forces along the bone, the normal and tangential components of the joint 
reaction forces, and the moments at the joints. This data is accessible as output 
of the dynamic gait model. 



CHAPTER VIII 
CONCLUSIONS 


The goal of this dissertation was to define a new method by which the ef- 
fects of musculature on stress development in the long bones can be ascertained 
throughout a normal human activity such as walking. As explained in Chapter II, 
it was a contention of this research that a forward analysis would provide a means 
by which the instantaneous loading conditions, such as joint reaction forces, joint 
moments, and muscular forces, of bone could be derived throughout the gait cy- 
cle. Since the musculotendon units are central to this approach, it was imperative 
that a good muscle model be utilized, one that reflected the basic characteristics 
of muscle while maintaining enough simplicity to render it useful in a simulation. 
Chapter III explained the key elements that a muscle model should incorporate 
while Chapter IV presented the model we felt was appropriate for motion simu- 
lations. Once the muscle model was determined, we sought a representation of 
the human body which included enough complexity to mimic the central issues of 
human gait. It was thought that the model developed by Yamaguchi succeeded 
in this regard. Chapter V illustrated the construction the complete model and 
the incorporation of the muscles as actuators in the system. Then Chapter VI 
reflected on the performance of the model in simulating gait. We felt that the 
simulations were indeed representative of human gait and that the derived loading 
conditions were in a reasonable range. Thus we proceeded in Chapter VII to do a 
finite element analysis of the tibia based on these results. Although modeling the 
bone as a uniaxial rod is quite simplified and only applicable during certain phases 
of gait, it still succeeded in showing the relevance of the musculature on the stress 
development in bones. 

In conclusion, a method for including the effects of muscular forces in a contin- 
uum analysis of the long bones has been presented and illustrated for a simulation 
of gait. The approach is multi-disciplinary in nature in that it hopes to bring 
together the often separately researched fields of muscle mechanics, bone mechan- 
ics and motion analysis. Although the example presented in this dissertation is 
illustrative of the method, there are several aspects of the research which should 
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be improved before the model will have any real predictive value in accessing the 
failure properties in bone. 

A list of the major improvements which would complement this dissertation 
follow: 

• A complete validation of the model needs to be performed. It is known 
that EMG activities supports the use of the muscles utilized in this analysis, 
and that the joint torques are in qualitative agreement with other torques 
presented in the literature. However, issues such as the trajectory of the 
center of mass and the joint reaction forces should be analyzed. 

• A more complex model of the body could be constructed. It would be advan- 
tageous to have a model which is capable of simulating the whole gait cycle. 
Another important improvement would be to break the torso segment into 
two segments so that the tilting of the pelvis could be incorporated into gait. 
There is a long list of possible complexities which could be added to make 
the simulated gait more representative of human gait; however, a balance 
should be maintained between the advantages of these complexities and the 
computational cost associated with their inclusion. 

• Explore the implications of more complex and realistic models of musculo- 
tendon dynamics. In particular the effects of muscle mass, muscle viscosity, 
and contractile series elasticity should be examined. Modulation of these 
components by activation, especially the contractile series element, should 
be investigated. 

• Incorporate additional muscle groups into the gait model. A major obstacle 
in this regard will be the computation of the neurological inputs required 
to control the model. It is anticipated that optimization techniques will 
be essential to this effort. A computationally efficient method called the 
pseudo-inverse method (Yamaguchi et al., 1995) has recently been developed 
for performing dynamic optimizations of movement. This method should be 
explored. 
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• Although the above improvements would be beneficial to the analysis, the 
most critical improvement will be in the implementation of a more complex 
model of the bone. One possible direction would be to follow the work of 
Toridis or Rybicki and modeling the bone as a beam which is transversely 
isotropic. This would allow for the inclusions of shear deformations and 
bending. It is also a model appropriate throughout the gait cycle, so that full 
use could be made of the derived loading conditions. Another direction would 
be to treat the whole skeletal structure as a space frame. Whatever model 
is utilized a careful interpretation of the output of the dynamic gait model 
is required to correctly interpret the loading conditions that are essential in 
the boundary value problems. 

• Characterize failure or damage in the skeletal system in terms of loading con- 
ditions associated with muscular and joint inputs. Two important issues that 
should be emphasized Eire muscle fatigue and pathological innervation. The 
ultimate objective of this work would be to relate critical stress to muscular 
failure. 
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% File 
% Problem 
% Desription 
% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% Note 


walk.al [AUTOLEV 3] 

7 link, 8 DOF model 

This model will have torques at all joints due to tendons, 
and ground reaction forces on both feet. It has two 
auxiliary coordinates to produce reaction forces at the 
shank joints in the normal direction. 

Stance-side muscles added: Soleus, Gastrocnemius, Vasti, 
Iliopsoas, and the Gluteus Medius/Minimus. 

Swing-side muscles added: Dorsiflexors, Vasti, 

Iliopsoas, and the Gluteus Medius/Minimus. 

Tendon and activation dynamics are included. 

All muscle functions will be added to the C-code as 
functions. Thus all changes to code can be applied 
directly. 

All muscle values of Yamaguchi are used explicitly. 

For muscle crossing the knee, the moment arm integration 
technique is used to figure Lmt, Vasti, Ham, Gastro. Since 
the Gastro is biarticulate across the knee and ankle, I used 
vector addition to find Lmt, but I did utitlize the moment 
arm data across the knee. 

All constants are entered. 


% 

% Physical Declarations: Newtonian, Frames, Particles, Bodies, Points 


% 

autoz on 
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newtonian N 

% Newtonian reference frame 

frames dd,ssft,sst,ssf,ssp 

% intermediate frames 

frames wwft,wwt,wwf,wwp 

% intermediate frames 

bodies a,b,c,d,e,f,g 

% 7 link model 

frames sft,st,sf,sp 

% reference frame for stance muscles 

frames wft,wt,wf,wp 

% reference frame for swing muscles 

points ab,ba,bc,cb,cd,dc 

% adjoining points 

points de,ed,ef,fe,fg,gf 

% adjoining points 

points ha,hg 

% heels of stance and swing foot 

points ta,tg 

% toes of stance and swing foot 

% 



% Muscle points in the stance leg 


% 


points so, si 
points gao,gai 
points vso,vsi,vsei 
points iso,isi,iseo 

points gmso,gmsi 
% 

% soleus 
% gastrocnemius 
% vasti 
% iliopsoas 

% gluteus medius/minimus 

/ v * — — — 

% Muscle points in swing leg 
% 

/ u 

points doo,doi,doeo 
points ho, hi 
points vwo,vwi,vwei 
points iwo,iwi,iweo 
points gmwo,gmwi 
% 

% dorsiflexors 
% hamstrings 
% vasti 
% iliopsoas 

% gluteus medius/minimus 

SU 

% Mathematical Declarations: 
% 

Variables, Constants, Specified, Mass 

/yj 

variables ul2’ 
variables q8’ 
variables FS\FGA’,FVS’ 

% generalized speeds(8) 

% generalized coordinates; derivatives 
% Force of muscles in stance leg 



84 


variables FIS’.FGMS’ 
variables FDO\FVW\FH’ 
variables FIW’,FGMW’ 
variables lmth’,lmtvs’,lintvw’ 
constants g 
constants w 

constants Ial,la3,llal,lla3 
constants lb, lib 
constants lc,llc 
constants ld.lld 
constants le=lc,lle=lc-llc 
constants lf=lb,llf=lb-llb 
constants Igl=la3,lg3=lal 
constants llgl=lla3, llg3=lal-llal 
mass a=ma,b=mb,c=mc,d=md,e 
inertia a, ial,ia2,ia3 
inertia b, ibl,ib2,ib3 
inertia c, icl,ic2,ic3 
inertia d, idl,id2,id3 
inertia e, iel=icl,ie2=ic2,ie3=ic3 
inertia f, ifl=ibl,if2=ib2,if3=ib3 
inertia g, igl=ial,ig2=ia2,ig3=ia3 

% 

% Geometry Relating the Segments 

% 

simprot(n,a,-2,ql) 

simprot(n,b,-2,q2) 

simprot(n,c,-2,q3) 

simprot(n,dd, 1 ,q4) 

simprot(dd,d,-2,q5) 

simprot(dd,e,2,q6) 

simprot(dd,f,2,q7) 


% Force of muscles in swing leg 

% Lengths of muscles spanning the knee 
% gravity 
% total weight 
% stance-leg foot 
% stance-leg shank 
% stance-leg thigh 
% trunk 

% swing-leg thigh 
% swing-leg shank 
% swing-leg foot 

me,f=mf, g=mg 
% stance-leg foot 
% stance-leg shank 
% stance-leg thigh 
% trunk 

% swing-leg thigh 
% swing-leg shank 
% swing-leg foot 


% pelvic list 
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simprot(dd,g,2,q8) 

simprot(a,ssft,2,pi-34*pi/180) % stance foot 

simprot(ssft,sft,l,pi/2) 
simprot(b,sst,-3,pi/2) 
simprot(sst,st,2,pi/2) 
simprot(c,ssf,-3,pi/2) 
simprot (ssf,sf ,2 ,pi/2 ) 
simprot (d , ssp ,-3 ,pi/2 ) 
simprot(ssp,sp,2,pi/2) 
simprot(d,wwp,-3,pi/2) 
simprot(wwp,wp,2,pi/2) 
simprot (e , w wf , 3 , pi/ 2 ) 
simprot (wwf,wf, -2 ,pi/2) 
simprot(f,wwt,3,pi/2) 
simprot(wwt,wt,-2,pi/2) 
simprot(g,wwft,-2,34*pi/180+pi/2) % swing foot 

simprot(wwft,wft,l,pi/2) 

% 

% Set up geometry for stance foot ”a” 

% 

p_ta_ab> =lal *al > 
p_ab_ha>=la3*a3> 

p_ta-ao>=llal*al>+lla3*a3> % center of mass 

% 

% Set Up Geometry for Stance Shank ”b” 

% 

p_ab.ba>=0> 
p_ba_bc>=lb*bl > 

p_ba_bo>=llb*bl> % center of mass 

% 

% Set Up Geometry for Stance Thigh ”c” 

% 


% stance tibial 
% stance femoral 
% stance pelvic 
% swing pelvis 
% swing femoral 
% swing tibial 
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p_bc_cb>=0> 

p_cb_cd>=lc*cl> 

p_cb_co>=llc*cl> % center of mass 

% 

% Set Up Geometry for Stance Trunk ”d” 

% 

p_cd_dc>=0> 
p_dc.de> =ld*d2> 

p_dc_do>=.5*ld*d2>+lld*dl> % center of mass 

% 

% Set Up Geometry for Swing Thigh ”e” 

% 

p_de_ed>=0> 

p_ed_ef>=le*el> 

p_ed_eo>=lle*el> % center of mass 

% 

% Set Up Geometry for Swing Shank ”f’ 

% 

p_ef_fe>=0> 
p _fe_fg> =lf*fl > 

pJe_fo>=llf*fl> % center of mass 

% 

% Set up Geometry for Swing Foot ”g” 

% 

p_fg_gf>=0> 

P-gf-hg>=lgl*gl> 

p_gf_tg>=lg3*g3> 

P-gf^go>=llgl*gl>+llg3*g3> % center of mass 

% 

% Soleus Stance Muscle 

% 

p_ba_so> =-.0292*st 1 > + .2467*st2> +.0006*st3> 



p_ba_si>=-.0365*sftl>-.0428*sft2>+.0056*sft3> 

% 

% Gastrocnemius Stance Muscle 

% 

p_cb-gao>=-.0203*sfl>+.0071*sf2>-.0073*sf3> 

p_ba_gai>=-.0368*sftl>-.0429*sft2>+.0028*sft3> 

% — 

% Vasti Stance Muscle 

% 

p_cb.vso> = .0 106*sf 1 > + .2026*sf2 > + .0205 *sf3 > 
p-ba_vsi>=.0170*stl>+.3930*st2>-.0006*st3> 

% 

% Iliopsoas Stance Muscle 

% 

p_cdJso>=.0075*spl>+.1350*sp2>-.04*sp3> 

P _bcJsi>=-.0180*sfl>+.3351*sf2>+.0116*sf3> 

p_cdJseo>=.0260*spl>+.0293*sp2>-.0042*sp3> 

% 

% Gluteus Medius Stance Muscle 

% 

p_dc_gmso>=-.0155*spl>+.0785*sp2>+.0076*sp3> 

p_cb^msi>=-.0159*sfl>+.3873*sf2>+.0589*sf3> 

% 

% Dorsiflexors Swing Muscle 

% — 

P-fg_doo>=-.0155*wtl>+.2175*wt2>+.0134*wt3> 

p_fg_doi>=.1035*wftl>-.052*wft2> 

p_fg_doeo>=.0259*wtl>+.0117*wt2>-.0093*wt3> 

% — 

% Hamstring Swing Muscle 

% 

p_de_ho>=-.0409*wpl>-.0455*wp2>-.014*wp3> 
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pjg-hi > =- . 0 1 7* wt 1 > + • 38 * wt 2 > + .0073 * wt3 > 

% 

% Vasti Swing Muscle 

% 

p_ef_vwo> =.0106* wfl > + . 2026*wf2> + .0205*wf3> 
p_fg_vwi>=.0170*wtl>-f .3930*wt2>-.0006*wt3> 

% 

% Iliopsoas Swing Muscle 

% 

p_de Jwo>= .0075*wpl > + . 1350*wp2 >- .04*wp3> 
p.ef Jwi>=-.0180* wfl > + .3351 * wf2> + .01 16*wf3> 
p_deJweo>=.0260*wpl>+.0293*wp2>-.0042*wp3> 

% 

% Gluteus Medius Swing Muscle 

% 

p_de_gmwo>=-.0155*wpl>+.0785*wp2>+.0076*wp3> 

p_ef-gmwi>=-.0159*wfl>+.3873*wf2>+.0589*wf3> 

% 

% Kinematic Differential Equations 

% 

ql’=ul 

q2’=u2 

q3’=u3 

q4’=u4 

q5’=u5 

q6’=u6 

q7’=u7 

q8’=u8 

% : 

% Angular Velocity of the Segments 

% 

w_a_n>=-ql’*n2> 
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w_b_n>=-q2’*n2> 
w_c_n > =-q3 ’ * n2 > 
w_dd_n>=q4’*nl> 
w_d_dd > =-q5 ’ * dd2 > 
w_d_n> =w_d_n> 
w_e_dd> =q6’*dd2 > 
w_e_n> =w_e_n> 
w_Ldd>=q7’*dd2> 
w_f_n > = w_f_n > 
w_g_dd>=q8’*dd2> 
w_g_n > = w-g_n > 

% 

% Velocities of Points 

% 

vta_n>=0> 

v2pts(n,a,ta,ab) 

v2pts(n,a,ta,ao) 

v2pts(n,a,ab,ha) 

% 

v_ba_n > =■ v_ab_n > 

v2pts(n,b,ba,bo) 

v2pts(n,b,ba,bc) 

% 

v_cb_n>=v_bc_n> 

v2pts(n,c,cb,co) 

v2pts(n,c,cb,cd) 

% 

v_dc_n > = v_cd _n > 

v2pts(n,d,dc,do) 

v2pts(n,d,dc,de) 

% 

v_ed_n>=v_de_n> 
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v2pts(n,e,ed,eo) 

v2pts(n,e,ed,ef) 

% 

vJe_n>=v_ef_n>+u9*fl>+ ul0*f3> 

auxiliary [1] =u9 

auxiliary[2]=ul0 

v2pts(n,f,fe,fo) 

v2pts(n,f,fe,fg) 

% 

v_gf_n>=vJg_n>+ull*fl>+ul2*f3> 

auxiliary [3] =u 1 1 

auxiliary [4] = u 1 2 

v2pts(n,g,gf,go) 

v2pts(n,g,gf,tg) 

v2pts(n,g,gf,hg) 

% 

% Motion Constraints 

% 

constrain ( auxiliary [u9 , u 1 0 , u 1 1 , u 1 2] ) 

% 

% Forces 

% 

Gravity(-g*n3>) 

% 

% Ground Reaction Forces and Torques at Heels 

% 

specified zheela’,zheelg’, delta’ 

% The following functions will needed to be added to the C-code. 
specified fa3,fgl,fg3,tt 
variables tabt,tbct,tcdt,tdet,teft,tfgt 
Variables rfe3,rfg3 
% 
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% The following functions will not be embedded in Z variables. 
zee_not=[tabt,tbct,tcdt,tdet,teft,tfgt,fa3,fgl,fg3,tt,rfel,rfe3,rfgl, 
rfgS.FS^GA.FVS^IS.FGMS^GMW.FIW^VW.FH.FDO] 

% 

zheela=dot (p_ta_ha> ,n3 > ) 
zheela’ =dt (zheela) 
zheelg=dot(p_ta_hg> ,n3>) 
zheelg’=dt(zheelg) 
heelnl=dot(v _hg_n> ,nl >) 
delta=-q8+2.164 
delta’=dt (delta) 

% 

% Ground reaction force on stance toe. 

Force_ha>+=fa3*n3> 

% 

% Ground reaction force on swing heel. 

Force_hg > + =fg 1 * n 1 > + fg3 * n3 > 

% 

% Torque on swing foot. 

Torque_g>+=tt*dd2> 

% 

% Passive Joint Moments 

% 

% The following variables specify the passive joint moments. 
Constants kab4,kbc4,kcd4,kde4,kef4 ) kfg5 
Constants theta_ab2, theta_bc2 ,theta_cd2 
Constants theta_de2,theta_ef2 ) theta_fg2 
Constants cab 1 ,cbc 1 ,ccd 1 ,cde 1 , cefl ,cfg 1 
% 

% Joint singles defined. 
theta_ab=ql-q2 +(34*pi/ 180) 
theta_bc=q2-q3 



theta_cd=q3-q5 

theta_de=pi-(q5+q6) 

theta_ef=q6-q7 

theta_fg=(124*pi/180)-(q8-q7) 

% 

% Passive torques defined. 

tabt=(kabl*exp(-kab2*(theta^b-theta_ab2))- 

kab3*exp(-kab4*(theta_abl-theta^ab))-cabl*dt(theta_ab)) 

tbct=(kbcl*exp(-kbc2*(theta_bc-theta_bc2))- 

kbc3*exp(-kbc4*(theta_bcl-theta_bc))-cbcl*dt(theta_bc)) 

tcdt=(kcdl*exp(-kcd2*(thetajcd-theta_cd2))- 

kcd3*exp(-kcd4*(theta_cdl-theta_cd))-ccdl*dt(theta.cd)) 

tdet=(kdel*exp(-kde2*(theta_de-theta_de2))- 

kde3*exp(-kde4* ( theta.de l-theta.de) )-cdel *dt ( theta.de) ) 

teft=(kefl*exp(-kef2*(thetajef-theta_ef2))- 

kef3*exp(-kef4*(thetajefl-theta_ef))-cefl*dt(theta.ef)) 

tfgt=(kfgl *exp(-kfg2* (theta Jg-thetaJg2) )- 

kfg3*exp(-kfg4* ( theta Jg 1- theta Jg) ) -cfg 1 * dt ( theta Jg) ) 

% 

torque(a/b,tabt*n2>) 
torque(b/c, tbct * n2 > ) 
torque.c>-=tcdt*n2> 
torque_d>+=tcdt*dd2> 
torque(e /d, tdet *dd2 > ) 
torque(f/e,teft*dd2>) 
torque(g /f, tfgt *dd2 > ) 

% 

% 

% Muscles 

% 

% 

% None of the lengths nor forces are normalized. 
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% Constants are input. 

% These functions will be specified in the C code 

% 

Specified fla,flp,fv,K,cs,Lt,Lm,mrv,mrh,mrg 
% time constant 
Constants tau 
% 

% Soleus muscle on stance leg 

% 

Specified lmtso’jlmsojltso.lwsojvmtsojvmsojvtsojcssojaso.faso.ms 
Constants F0mso,L0mso,Lstso,Paso 
% 

ms=dot(cross(p_ab-SO>,unitvec(p_so-si>)),n2>) 
lwso= LOmso*sin(Paso) 
lmtso=mag(p-SOJ3i>) 
vmtso=dt(lmtso) 
ltso=Lstso*Lt 
lmso=Lm 
csso=cs 
% 

faso=( (FS/csso) - F0mso*flp) / ( FOmso*aso*fla ) 
vmso= (LOmso/tau)*fv/csso 
vtso= vmtso- vmso 
FS’=(FOmso/lstso) :,! K*(vtso) 

% 

torque(a/b,FS*ms*n2>) 

% 

% Gastrocnemis muscle on stance leg 

% 

Specified lmtga’,lmga,ltga,lwga,vmtga,vmga,vtga,csga,aga,faga,mgak,mgaa 
Constants FOmga,LOmga,Lstga,Paga 
% 


mgak=-mrg*(q3-q2) 

mgaa=dot(cross(p^ab.gao> ,unitvec(p_gao^ai> )) ,n2>) 

1 wga=LOmga* sin( Paga) 
lmtga=mag(p-gao.gai>) 
vmtga=dt (lmtga) 
ltga=Lstga*Lt 
lmga=Lm 
csga=cs 
% 

faga=((FGA/csga)-FOmga*flp)/(FOmga*aga*fla) 
vmga=(LOmga/tau)*fv/csga 
vtga=vmtga-vmga 
FGA’=(FOmga/Lstga) *K* ( vtga) 

% 

torque(a/b,FGA*mgaa*n2>) 

torque(b/c,FGA*mgak*n2>) 

% 

% Vasti muscle on stance leg 

% 

Specified lmvsjltvs^wvs.vmtvsjvmvs^tvs^svs^vs^avs.mvs 
Constants FOmvs,LOmvs,Lstvs,Pavs 
% 

mvs=mrv* (q3-q2) 
lwvs=LOmvs*sin(Pavs) 

Imtvs’==inrv*(q3-q2)*(u3-u2) 

vmtvs=lmtvs’ 

ltvs=Lstvs*Lt*(lmtvs-lmvs) 

lmvs=Lm 

csvs=cs 

% 

favs=((FVS/csvs)-FOmvs*flp) /(FOmvs*avs*fla) 
vmvs=(LOmvs/tau)*fv/csvs 


vtvs=vmtvs-vmvs 

FVS’=(FOmvs/Lstvs)*K*(vtvs) 

% 

torque(b/c,FVS*mvs*n2>) 

% 

% Iliopsoas muscle on stance leg 

% 

Specified lmtis , ,lmis,ltis,lwis,vmtis,vmis,vtis,csis,ais,fads 
Specified misddl,misdd2,misn2 
Constants FOmis,LOmis,Lstis,Pais 
% 

misdd2=dot (cross(p_cd Jseo> ,unitvec(p Jseo Jsi> )) ,dd2 >) 
misddl=dot(cross(p_cd_iseo>,unitvec(piseoJsi>)),ddl>) 
misn2=dot (cross(p.cd Jseo> ,unitvec(p Jseo Jsi >) ) ,n2>) 
lwis=LOmis*sin(Pais) 

lmtis=(mag(pJsoJseo>)+mag(p_iseoJsi>)) 
vmtis— dt(lmtis) 
ltis=Lstis*Lt 
lmis=Lm 
csis=cs 
% 

fais= ((FIS / csis)-FOmis*flp) /(FOmis* AIS *fla) 
vmis=(LOmis/tau)*fv/csis 
vtis=vmtis-vmis 
FIS’=(F0mis/lstis)*K*(vtis) 

% 

torque_c> + =-FIS *misn2*n2 > 
torque_d>+=FIS*misdd2*dd2> 
torque_d> + =FIS * misdd 1 * dd 1 > 

% 

% Gluteus Medius/Minimus muscle on stance leg 

% 


Specified lmtgms’,lmgms,ltgms,lwgms,vmtgms,vmgnis,vtgms 
Specified csgms,agms,fagms 
Specified mgmsdd2,mgmsddl,mgmsn2 
Constants FOmgms,LOmgms,Lstgms,Pagms 
% 

mgmsdd2=dot (cross ( p.cd _gmso >, uni tvec (p-gmso-gmsi >)) ,dd2 > ) 

mgmsddl=dot(cross(pjcd_gmso>,unitvec(p_gmso_ginsi>)),ddl>) 

mgmsn2=dot(cross(pxd_gmso>,unitvec(p_gmso-gmsi>)),n2>) 

lwgms=LOmgms*sin(Pagms) 

lmtgms=mag(p_gmso-gmsi> ) 

vmtgms=dt (lmtgms) 

ltgms=Lstgms*Lt 

lmgms=Lm 

csgms=cs 

% 

fagms=((FGMS/csgms)-FOmgms*flp) / (FOmgms*AGMS*fla) 
vmgms= (LOmgms / tau) *fv/ csgms 
vtgms=vmtgms-vmgms 
FGMS’=(FOmgms/lstgms)*K*(vtgms) 

% 

torque_c>+ =-FGMS *mgmsn2 *n2 > 
torque_d>+=FGMS*mgmsdd2*dd2> 
torque.d > += FGMS * mgmsdd 1 * dd 1 > 

% 

% Dorsiflexors muscles on swing leg 

% 

Specified lmtdo’,lmdo,ltdo,lwdo,vmtdo,vmdo,vtdo,csdo,ado,fado,mdo,fsdo 
Constants F0mdo,L0mdo,Lstdo,Pado 
% 

mdo=dot(cross(p_fg_doeo>,unitvec(p_doeo_doi>)),dd2>) 
lwdo=L0mdo* sin( Pado) 

lmtdo=(mag(p-doo.doeo>)+mag(p_doeo-doi>)) 


97 


vmtdo=dt(lmtdo) 

ltdo=Lstdo*Lt 

lmdo=Lm 

csdo=cs 

% 

fado=((FDO/csdo)-FOmdo*flp)/(FOmdo*ADO*fla) 
vmdo= ( LOmdo / tau) *fv/csdo 
vtdo=vmtdo-vmdo 
FDO’=(FOmdo/Lstdo)*K*(vtdo) 

% 

fsdo=(FDO)*dot(fl>,unitvec(p-doeo_doi>)) 

% 

torque(g/f,FDO*mdo*dd2>) 

% 

% Hamstring muscle in swing leg 

% 

Specified lmh.lthdwh.vmt^vm^vthjCsh^h.fa^mhk.mhh.fsh 
Constants F0mh,L0mh,Lsth,Pah 
% 

mhk=mrh*(q7-q6) 

mhh=dot(cross(p_de_ho>,unitvec(p_ho_hi>)),dd2>) 

lwh=LOnih*sin(Pah) 

Imth’=mrh*(q7-q6)*(u7-u6) 

vmth=lmth’ 

lth=Lsth*Lt*(lmth-lm) 

lmh=Lm 

csh=cs 

% 

fah=((FH/csh)-FOmh*flp)/(FOmh*AH*fla) 

vmh=(LOmh/tau)*fv/csh 

vth=vmth-vmh 

FH’=(FOmh/Lsth)*K*(vth) 



% 

fsh=FH*dot(fl>,unitvec(p_ho_hi>)) 

% 

torque(f/e,FH*mhk*dd2>) 

torque(e/d,FH*mhli*dd2>) 

% 

% Vasti muscle on swing leg 

% 

Specified lmvw,ltvw,lwvw,vmtvw,vmvw,vtvw,csvw,avw,favw,mvw,fsvw 
Constants FOmvw,LOmvw,Lstvw,Pavw 
% 

mvw=mrv* (q7-q6) 
lwvw=LOmvw*sin (Pavw) 

Imtvw’=mrv*(q7-q6)*(u7-u6) 

vmtvw=lmtvw’ 

ltvw=Lstvw*Lt*(lmtvw-lmvw) 

lmvw=Lm 

csvw=cs 

% 

favw=((FVW/csvw)-FOmvw*flp) / (FOmvw*avw*fla) 
vmvw=(LOmvw/tau) *fv/csvw 
vtvw=vmtvw-vmvw 
FV r W’=(FOmvw/Lstvw)*K*(vtvw) 

% 

fsvw=FVW*dot (fl > ,unitvec(p_vwo_vwi>)) 

% 

torque(f/e,FVW*mvw*dd2>) 

% 

% Iliopsoas muscle on swing leg 

% 

Specified lmtiw , ,lmiw,ltiw,lwiw,vmtiw ) vmiw,vtiw,csiw,aiw > faiw,miw 
Constants FOmiw,LOmiw,Lstiw,Paiw 
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% 

miw=dot(cross(p_de Jweo> ,unitvec(p Jweo Jwi> ) ) ,dd2>) 
lwiw=LOmiw*sin(Paiw) 

lmtiw= ( mag( pJ wo Jweo >)+ mag(p Jweo J wi > ) ) 

vmtiw— dt ( lmti w) 

ltiw=Lstiw*Lt 

lmiw=Lm 

csiw=cs 

% 

faiw=((FIW/csiw)-FOrmw*flp)/(FOiniw*AIW*fla) 
vmiw= (LOmi w/tau) *fv/ csiw 
vtiw=vmtiw-vmiw 
FIW’=(FOmiw/Lstiw)*K*(vtiw) 

% 

torque(e/d,FIW*miw*dd2>) 

% 

% Gluteus Medius/Minimus muscle on swing leg 

% 

Specified lmtgmw’.lmgmwjltgmwjlwgmwjvmtgmwjvmgmw^gmw 
Specified csgmw,agmw,fagmw,mgmw 
Constants FOmgmw,LOmgmw,Lstgmw,Pagmw 
% 

mgmw=dot(cross(p_de-gmwo>,unitvec(p_gmwo.gmwi>)),dd2>) 

% 

lwgmw=LOmgmw*sin(Pagmw) 

lmtgmw= mag(p-gmwo _gmwi > ) 

vmtgmw=dt (lmtgmw) 

ltgmw=Lstgmw*Lt 

lmgmw=Lm 

csgmw=cs 

% 

fagmw=((FGMW/csgmw)-FOmgmw*flp)/(FOmgmw*AGMW*fla) 


100 


vmgmw = (LOmgmw / tau) *fv / csgmw 

vtgmw=vmtgmw-vmgmw 

FGMW’=(FOmgmw/Lstgmw)*K*(vtgmw) 

% 

torque(e/d,FGMW*mgmw*dd2>) 

% 

% Forces at the ankle and knee of the swing side shank 

% 

Force(ef/fe,rfel*fl>+rfe2*f2>+rfe3*f3>) 

Force(fg/gf,rfgl*fl>+rfg2*f2>+rfg3*f3>) 

% 

Contributions=fr() 

%__ “ 

% Segmental Torques 

% — — 

torque_a> 

torque_b> 

torque_c> 

torque_d> 

torque_e> 

torqueJF> 

torque_g> 

% 

% Equations of Motion 

% 

zero= fr()+frstar() 
kane(rfel ,rfe3,rfgl,rfg3) 

% 

% Expression to be Output by the C Code 


% 

output T,ql 
output T,q2 
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output T,q3 
output T,q4 
output T,q5 
output T,q6 
output T,q7 
output T,q8 
output etc. 

% 

% Units Constants for C Code 

% 

units [t]=s 
units [tau]=s 

units [ql,q2,q3,q4,q5,q6,q7,q8]=deg 
units [ul,u2,u3,u4,u5,u6,u7,u8]=r/s 

units [FS,FGA,FVS ) FIS,FGMS > FDO,FVW,FH,FIW,FGMW]=N 
units [lal,la3,Ual ,llci3,lb,llb,lc,llc ) ld,lld]=m 
units [ma ) mb,mc,md,me,mf,mg,w]=kg 
units [ial,ia2,ia3,ibl,ib2,ib3,icl,ic2,ic3,idl,id2,id3]=kg*tn2 
units [cabl,cbcl,ccdl,cdel,cefl,cfgl]=(N*s)/(r*m) 

units [kabl,kab3,kbcl,kbc3,kcdl,kcd3,kdel,kde3,kefl,kef3,kfgl,kfg3]=N*m 
units [kab2,kab4,kbc2,kbc4,kcd2,kcd4,kde2 ) kde4,kef2,kef4,kfg2,kfg4]=l/r 
units [theta-abl,theta_ab2,theta_bcl,theta_bc2,theta_cdl,theta_cd2]=r 
units [theta.del, theta_de2,theta_efl,theta_ef2,theta_fgl,theta_fg2]=r 
units [g]=m/s2 

units [F0mSo,F0mGA,F0mVS,F0mIS,F0mGMS]=N 
units [FOmDO,FOmVW,FOmH,FOmIW,FOmGMW]=N 
units [L0mSo,L0mGA,L0mVS,L0mIS,L0mGMS]=m 
units [LOmDO,LOmVW,LOmH ) LOmTW ) LOniGMW]=m 
units [LstSo,LstGA,LstVS,LstIS,LstGMS]=m 
units [LstDO,LstVW,LstH,LstIW,LstGMW]=m 
units [PaSo,PaG A, PaVS ,PaIS ,PaGMS] =deg 
units [PaDO ,PaVW,PaH,PaIW,PaGMW] =deg 



input g=9.81 

input ial=.002,ia2=.008,ia3=.009,ibl=.019,ib2=.065,ib3=.065 
input icl=.080,ic2=.126,ic3=.126,idl=.764,id2=3.407,id3=3.297 
input lal=. 175, la3=. 118, lb=. 435, lc=.410,ld=. 172, llal=.l,Ua3=. 0295 
input Ub=.247,Dc=.227,lld=.343,ma=l.l,mb=3.75,mc=7.58,md=5L22 
input me=7.58,mf=3.75,mg=l.l 

input cabl=.943,cbcl=3.17,ccdl=1.09,cdel=1.09,cefl=3.17,cfgl=.943 

input kabl=2,kab2=5,kab3=9,kab4==5 

input kbcl=3.1,kbc2=5.9,kbc3=10.5,kbc4=11.8 

input kcdl=2.6,kcd2=5.8,kcd3=8.7,kcd4=1.3 

input kdel=2.6,kde2=5.8,kde3=8.7,kde4=1.3 

input kefl=3.1,kef2=5.9,kef3=10.5,kef4=11.8 

input kfgl=2,kfg2=5,kfg3=9,kfg4=5 

input theta_ab 1=1.92, theta_ab2 =1.047, theta_bc 1 =0 , theta_bc2=- 1.92 

input theta_cdl=1.92,theta_cd2=.1744,theta_del=1.92,theta_de2=.1744 

input theta_efl=0,thetajef2=-1.92,theta_fgl=1.92,thetaJg2=1.047 

input F0mso=3599,L0mso=.0243,Lstso=.27,Paso=25 

input F0mga=1423, L0mga=. 0482, Lstga=. 4250, Paga= 14.8 

input F0mvs=6482,L0mvs=. 1096,Lstvs=.2250,Pavs=4.5 

input F0mis=1474,L0mis=. 1269, Lstis=. 0850, Pais=7 

input F0mgms=2686,L0mgms=.0760,Lstgms=.0355,Pagms=10.4 

input F0mdo= 1400, L0mdo=. 1009, Lstdo=. 2250, Pado=6.9 

input F0mvw=6482,L0mvw=.1096,Lstvw=.2250,Pavw=4.5 

input F0mh=2348,L0mh=.1065,Lsth=.3850,Pah=8.7 

input F0miw=1474,L0miw=. 1269, Lstiw=. 0850, Paiw=7 

input F0mgmw=2686,L0mgmw=.0760,Lstgmw=.O355,Pagmw=10.4 

input tinitial=0,tfinal=.6,integstp=.01 

% 

% C Code generation for numerical solution 
digits 6 

code dynamics() walk.c, subs 


